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1.  Introduction 


The  determination  of  the  natural  modes  of  vibration  in  layered  elastic  media  is  a  central  problem 
of  geophysics  since  earthquake-induced  free  oscillations  of  the  Earth  are  often  used  to  constrain 
models  of  the  Earth’s  interior.  Poisson  (1)  was  the  first  to  determine  the  free  radial  vibrations  of 
a  homogeneous  sphere,  and  Lamb  (2)  calculated  some  of  its  natural  frequencies;  the  introduction 
in  Love’s  Mathematical  Treatise  of  Elasticity  (3)  provides  an  historical  account  of  the  early 
developments  in  this  field.  It  was  only  with  the  advent  of  high-speed  computers  that 
planetary-scale  vibration  problems  could  be  solved  with  finite,  time-dependent  sources  in 
compressible,  self-gravitating,  inhomogeneous,  and/or  layered  oblate  spheroids  of  revolution 
(4,  5).  Interestingly,  the  observed  resonance  in  the  free  oscillations  of  the  solid  Earth,  has  been 
attributed  to  fluid-structure  interaction  induced  by  the  atmosphere,  rather  than  the  cumulative 
effect  of  microearthquakes  (6). 

Studies  of  natural  vibrations  in  elastic  media  at  spatial  scales  of  technological  interest  focus  on 
resonance  in  layered  structures  (7),  anisotropic  elastic  bodies  (3),  anisotropic  layered  crystals 
(9,  10),  elastic  plates  (11,  12),  periodic  media  (13),  laminated  and  sandwich  plates  (14), 
composite  laminates  (15),  piezoelectric  composites  (16),  locally  resonant  acoustic  metamaterials 
(17,  18),  and  acoustic  wave  devices  that  exploit  the  phenomenon  (19-21).  Thus,  it  is  seen  that 
resonance  can  enhance  the  performance  of  many  sensors  and  devices,  yet  can  damage  structures 
subjected  to  sustained,  temporally  periodic  loading  (22). 

Although  there  are  numerous  works  that  establish  the  natural  modes  of  vibration  of  multilayered 
media  (23,  24),  and  analogous  problems  for  n-segmented  strings  (25)  and  n-strings  (26), 
analytical  solutions  for  the  transient  resonance  response  of  multilayered  media  have  been  limited 
to  analyses  involving  only  a  few  layers.  Churchill  (27)  uses  Laplace  transform  methods  to  derive 
the  transient  resonance  response  of  a  free-fixed  elastic  bar,  and  notes  the  presence  of  the  product 
of  a  temporal  term  and  time-harmonic  function  in  the  solution,  which  becomes  unbounded  at 
large  times;  similar  unbounded  behavior  is  observed  in  the  interaction  of  a  compressible  fluid 
with  an  elastic  plate  (28).  Caviglia  and  Morro  provide  closed-form,  time-domain  expressions  for 
transient  waves  in  a  multilayered  elastic  (possibly  anisotropic)  medium  (29).  They  derive  the 
transient  resonance  response  of  a  single  isotropic  elastic  layer  sandwiched  between  two 
half-spaces  that  is  subjected  to  a  temporally  periodic  sawtooth  function,  and  observe  that 
resonance  makes  the  elastic  layer  acoustically  transparent.  A  similar  two-dimensional  problem  is 
studied  by  Kaplunov  and  Krynkin  (30),  who  examine  the  influence  of  layer  stiffness  and 
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thickness  on  the  resonance  vibrations  in  a  symmetrically  point-loaded  single  layer  of  elastic 
material  enclosed  between  two  half-spaces.  Qiang  et  al.  (31)  study  the  natural  frequencies  of 
anisotropic  multilayers,  and  illustrate  beat  phenomena  in  two  and  three  layer  systems  using  an 
efficient  numerical  eigensolution  scheme  that  is  based  on  semi-analytical  methods.  Exact 
expressions  for  the  reflection  coefficients  for  a  two-dimensional  elastic  layer  overlying  an  elastic 
half-space  are  obtained  by  Fokina  and  Fokin  (32),  but  the  transient  response  of  the  system  at 
resonance  is  not  analyzed.  Graff  presents  a  general  method  (33)  to  determine  the  natural 
frequencies  of  composite  rods  for  power  ultrasonic  applications  with  a  specific  solution  for  a 
two-rod  system.  The  method  can  be  extended  to  multilayered  systems,  but  the  general  problem 
of  determining  the  natural  frequencies  with  specific  geometry  and  material  properties  can  only  be 
solved  using  numerical  methods.  Exact  analytical  expressions  are  derivable,  however,  for  the 
transient  resonance  response  of  a  multilayered  medium  if  we  specialize  the  medium  to  be  of 
Goupillaud-type  (34).  The  Goupillaud  specialization  is  often  used  in  geophysical  applications  to 
model  wave  propagation  in  inhomogeneous  media  (35).  The  present  treatment  uses  a  global 
matrix  method  that  is  attributed  to  Knopoff  (36),  rather  than  the  Thomsen-Haskell  transfer  matrix 
formalism  (37,  38).  Since  the  stresses  in  the  multilayered  medium  are  written  in  terms  of 
recursion  relations  that  interlace,  only  half  the  number  of  equations  are  required  relative  to  those 
in  classical  global  matrix  methods  (36).  Finally,  Bube  and  Burridge  (39)  treat  the  inverse 
problem  of  finding  the  coefficients  of  a  first  order  2x2  hyperbolic  system  related  to  reflection 
seismology.  In  order  to  numerically  solve  the  continuous  initial-boundary-value  problem,  several 
difference  schemes  are  applied  as  discretizations  to  the  corresponding  differential  equations.  The 
difference  scheme  IVp,  given  in  equation  (3.1.11),  pg.  5 17  of  Bube  and  Burridge  (39),  appears 
to  match  our  recursive  relations  for  the  stress  values  given  in  section  2.  Our  recursive  relations 
are  exact  for  discretely  layered  media  of  Goupillaud-type,  and  the  stress  solutions  that  are  derived 
herein  using  ^-transform  methods,  are  then  used  to  determine  the  resonance  response  for  the 
discrete  model,  and  subsequently  extended  to  the  continuous  model. 

The  remainder  of  the  report  is  organized  as  follows.  The  resonance  frequency  spectrum  for  the 
discrete  model  is  developed  in  section  2  by  generalizing  the  global  system  of  recursion  relations 
derived  in  Velo  and  Gazonas  (40).  Explicit  analytical  expressions  for  the  stress  in  a  multilayered 
elastic  medium  are  obtained  using  ^-transform  methods  and  are  quantitatively  compared  with  the 
recursion  relation  solutions  at  resonance.  These  results  are  then  extended  in  section  3  to 
analytically  describe  the  natural  frequency  spectrum  of  a  free-fixed  Goupillaud-type  strip,  also 
obtained  using  a  generalization  of  Graff’s  method.  Various  applications  and  interpretations  of  the 
frequency  results  are  discussed  in  section  4.  Section  5  summarizes  the  results,  while 
miscellaneous  references  and  derivations  are  included  in  the  appendices,  including  a  proof  that 
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the  zeros  of  the  palindromic  polynomial,  representing  the  determinant  of  the  global  system 
matrix,  lie  on  the  unit  circle  for  1  <  m  <  5. 


2.  Forced  Vibrations  and  Resonance  Frequencies:  The  Discrete  Model 


A  finite  strip  made  of  m, -layers  of  homogeneous  elastic  materials  is  considered,  where  the  density 
and  elastic  modulus  of  the  material  occupying  the  it h  layer  are  represented  by  p,  and  E, 
respectively,  for  i  —  1, . . . ,  m.  Using  the  definition  of  the  wave  speed  c  =  (E/p) l/2,  we  relate  to 
the  zth  layer  the  wave  speed  c,  for  ?'  =  1,2,...,  m.  The  characteristic  impedance  in  the  / 1 h  layer 
is  given  as  a  product  of  the  material  density  and  wave  speed  q,  while 
a,-  =  — — —  =  /EiPi  ,  represents  the  impedance  ratio  between  layers  i  and  /'  +  1  for 

Ci+lPi+l  yjEi+lPi  +  1  F  V  J 

i  —  1, . . .  ,m  —  1.  The  m-layered  strip  is  of  Goupillaud-type,  which  means  equal  wave  travel 
time  for  each  layer.  By  replacing  the  spatial  variable  x  of  the  wave  equation  with  the  new 
variable  £  =  [/  /j/-,  the  problem  can  be  converted  to  a  Goupillaud-type  strip  with  the  same  length 
L  as  the  transit  time  r  through  the  strip,  equal  layer  lengths  of  and  an  equal  wave  travel  time  of 
E  for  each  layer  in  either  direction,  which  further  implies  (equal)  wave  speed  of  unity  for  each 
layer.  Here  c  =  c(s )  is  the  piecewise  constant  wave  speed  function  taking  values  ci,  c2, . . . ,  cm  in 
each  layer,  respectively.  In  sections  2  and  3,  it  is  shown  that  the  stress  values  and  impedance 
ratios  remain  unchanged  as  a  result  of  this  simplification.  Therefore,  although  the  spatial  variable 
transformation  mentioned  previously  simplifies  the  problem,  and  is  used  to  analyze  a 
Goupillaud-type  strip  with  length  r  and  unit  wave  speed  in  each  layer,  the  results  can  be 
applied/interpreted  for  a  Goupillaud-type  strip  of  unequal  layer  lengths. 

In  prior  work  (40),  the  strip  was  subjected  to  zero  initial  conditions,  and  a  Heaviside  step  in  stress 
loading  p  =  constant  at  time  /  =  0  at  one  end,  and  held  fixed  at  the  other  end.  This  is  illustrated 
in  figure  1  with  f(n)  —  p  for  n  >  0.  The  present  work  generalizes  this  system  by  replacing  the 
Heaviside  step  in  stress  loading  with  a  discrete  loading  function  f(n),n  >  0.  After  transforming 
the  spatial  variable  x  of  the  wave  equation  to  the  new  variable  £  =  f*  /!/-.  the  initial-boundary 
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value  problem  becomes 


/  <92u(£,t)  _  d2u(^,t) 

dt2  ~  as,2 


<r(0,  t)  =  f|(0,i)  =  -Bi|(0,t) 


u(r,t)  =  0, 


for  ('1  <  £  <  — ,  i  =  1, . . . ,  m, 

m  ^  m 7 

f(n)  for  (n  —  1)—  <  t  <  n—,  n  >  1, 

J  \  /  \  j  m  m 7  —  ‘ 


a) 


l  «K,o)  =  S(5,o)  =  o. 

The  functions  u  and  a  represent  the  displacement  and  stress,  respectively,  while  f(n)  represents 
the  discrete  applied  loading.  Here,  Et  and  represent  the  actual  material  properties  for  the 
physical  case,  while  pi  and  Et  represent  the  material  properties  of  the  ?th  layer  for  the  simplified 
case.  We  will  use  this  notation  whenever  we  want  to  distinguish  the  physical  case  from  the 
simplified  case.  Having  p%  —  Et  —  —  —  \J E,pr,  it  implies  that  the  impedance  ratio 

_  EWE  _  \[Ei>i  _  _Ej_  as  wej]  as  t^e  stress  values  remain 

the  same  for  i  —  1,  2, . . . ,  m  —  1.  As  seen  in  figure  1,  the  previous  coordinates  of  the  boundary 
and  layer  interfaces  x0  =  0,  Xj,  and  xrn  =  L  are  replaced  by  the  new  coordinates  £o  =  0, 
and  =  r  for  j  =  1,2,....  rn.  The  time  variable  t  is  represented  on  the  vertical  axis  and  ^ 
represents  the  (equal)  wave  travel  time  for  each  layer  of  the  m-layered  strip  in  either  direction.  In 
figure  1,  the  stress  jump-discontinuities  are  propagated  along  the  dashed  characteristic  lines  of  the 
wave  equation,  while  the  inner  vertical  solid  lines  represent  the  layer  interfaces.  Due  to  the  zero 
initial  conditions,  the  stress  values  below  the  first  characteristic  line  segment  with  equation  t  =  £ 
for  0  <  £  <  r  are  zero.  Above  it,  however,  the  intersection  of  characteristics  with  each  other  and 
the  boundaries  splits  the  region  into  square  subregions.  Due  to  the  continuity  of  stress  and 
displacement  at  each  layer  interface,  the  stress  value  is  constant  in  each  square  subregion  and  is 
represented  by  Sj(n)  for  i  =  1,2,...,  rn  and  n  >  1.  In  a  given  layer  i,  the  stress  values  alternate 
(interlace)  between  Sj_i(n)  and  s*(n)  for  n  =  1,  2, . . .,  i  =  1,2, ...  ,m,  and  s0(n)  =  f{n).  When 
positioned  in  the  middle  of  the  ith  or  (i  +  l)th  layer,  the  time  interval  at  which  s^n)  is  reached  is 
4n  2m~3  <  t  <  in\ as  depicted  later  in  figures  2,  3,  4,  6,  8,  9,  and  11.  Here  n  >  1, 
i  =  0,1,2, ...  ,m,  and  s0(n)  =  f(n). 
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Figure  1 .  Lagrangian  diagram  for  an  elastic  strip  made  of  m-layers  of  equal  wave 
travel  time. 


From  the  superposition  of  the  transmitted  and  reflected  waves  at  the  layer  interfaces  and 
boundary,  a  system  of  coupled  first  order  difference  equations  for  the  stress  terms  was  derived  in 
Velo  and  Gazonas  (40);  the  system  was  then  solved  using  ^-transform  methods.  Herein,  we 
generalize  this  system  of  equations  by  incorporating  an  arbitrary  discrete  loading  function  /(n), 
n  >  0  which  appears  in  the  first  line  of  equation  2: 

'  si(n  +  1)  =  -si(n)  +  j^-s2(n)  +  ^/(n  +  1), 
s2(n  +  1)  =  -s2(n)  +  Yyt2s^(n)  +  j^si(n  +  1), 

<  Si{n  +  1)  =  - Si(n )  +  5^Si+i(ra)  +  ^s*_i(ra  +  1),  (2) 

sm-i(n  +  1)  =  — sm_i(n)  +  sm(n)  +  1+cfmi  sm_2(n  +  1), 

,  sm{n  +  1)  =  ~sm(n )  +  2sm-i(n  +  1), 
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subject  to  zero-stress  initial  conditions  Sj(0)  =  0  for  1  <  i  <  m.  Here  n  >  0  is  a  suitable  time 
index  and  Si(n)  represent  the  stress  terms  for  1  <  i  <  m  (see  figure  1).  The  discrete  function 
f(n),  n  >  0,  represents  the  stress  loading  applied  at  £  =  0.  Notice  that  by  applying  the  limiting 
condition  ct*  — *  0  to  the  recursive  relation  Si(n  +  1)  =  Si(n)  +  ^t.si+ i(n)  +  +  1) 

for  i  =  m,  we  recover  the  last  line  of  equation  2,  corresponding  to  the  fixed  end.  A  similar 
system  of  difference  equations  may  be  obtained  for  other  initial  and  boundary  conditions.  For 
example,  a  free-free  multilayered  system  can  be  studied  by  setting  the  last  line  in  equation  2  to 
sm(n)  =  sm( 0)  =  0,  and  a  multilayer  on  a  half-space  can  be  obtained  using 
sm(n  +  1)  =  sm-i(n  +  1),  and  the  rest  of  the  equations  for  1  <  i  <  m  remain  the  same. 

Equation  2  can  then  be  solved  using  ^-transform  methods.  Indeed,  due  to  the  zero-stress  initial 
conditions  Sj(0)  =  0  for  1  <  i  <  m,  the  z -transform1  of  S;(n  +  1)  becomes 

Z(si(n  +  1))  =  zSi(z)  -  zsi(Q)  =  zSi(z),  (3) 

where  St(z)  is  the  ^-transform  of  sr{n),  while 

Z(f(n  +  1))  =  zF(z)  -  zf( 0)  =  z(F(z)  -  /( 0)),  (4) 

where  F(z)  is  the  ^-transform  of  f(ri).  As  a  result,  the  ^-transform  of  equation  2  is  expressed  as 

'  zSt(z)  =  —Si(z)  +  ^S2(z)  +  ^z{F(z)  -  /( 0)), 

;  zS2(z)  =  —S2(z)  +  S3(z )  +  j^zS^z), 

<  zSi(z)  =  —Si(z)  +  ^:Si+1(z)  +  j^zSi-^z),  (5) 

zSm-i(z)  =  -Sm-i(z)  +  ^^Sm{z)  +  T+£^zSm-2{z), 

<  zSm^z')  =  SmiyZ')  T  2zSm-l(z). 


After  reorganizing  the  terms  in  equation  5,  the  linear  system  is  written  in  matrix- vector  form  as 


A  v  —  K 


(6) 


'The  definition  for  the  ^-transform  of  g(n),  n  >  0,  is  Z(g(n ))  =  G(z) 
plane  (see  Jury  (41)). 


OO 

E  9{n)z~n  for  |z|  >  R  in  the  complex 

n— 0 
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where  Am  is  a  tridiagonal  matrix  given  as 


Am 


z  +  1 

-77i«i 

0 

0 

0 

0 

0 

~V2Z 

z+1 

—  7/202 

0 

0 

0 

0 

0 

Z  +  1 

—r]3a3 

0 

0 

0 

0 

0 

0 

0 

~T]m-lZ  Z  +  1 

Vm—lC^m 

0 

0 

0 

0 

1 

o 

Z+1 

) 


mxm 


(7) 


while 


Si(z) 

'  mz(F(z)  -  /( 0))  ' 

S2(z) 

0 

,  bm  — 

Sm{z) 

0 

- 

mx  1 

. 

Vi  =  1^7  for*  =  l,rjm  =  2.  (8) 


Due  to  the  sparseness  of  bm,  the  solution  of  the  linear  system  in  equation  6  is 


(-1)1+1|Ai,!| 

1 

rH 

rH 

1 _ 

(— i)1+2  a1i2 

rH 

< 

mz(F(z)-f(0)) 

mz(F(z)-f(  0)) 

|Am| 

|Am| 

(— l)1+m  Aljm 

(-l)1+m|A1,m| 

(9) 


Here  |  Am|  is  the  determinant  of  Am,  Ai ,i  for  i  —  1, . . . ,  m  are  minors  of  Am,  and  |  Ax j|  are  their 
corresponding  determinants.  The  following  recursive  relation  (see  e.g.,  Muir  (42))  holds 


Am|  =  (z  +  1) | Am_i |  -  rjirj2aiz\Am_2l 


(10) 
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where  Am  corresponds  to  the  m- layered  strip  with  impedance  ratios  a  \ .  a2,  •  •  • ,  om- 1 ;  and  Am_! 
corresponds  to  the  remaining  (m  —  1) -layered  strip  with  impedance  ratios  a2,  a3, . . . ,  arn-  i , 
obtained  after  the  removal  of  the  first  layer  from  the  m-layered  strip;  Am_2  corresponds  to  the 
remaining  (m  —  2)-layered  strip  with  impedance  ratios  q3,  . . . ,  am_ , ,  obtained  after  the  removal 
of  the  first  two  layers  from  the  m-layered  strip.  One  can  derive  by  induction  that  the  determinant 
|  Am  |  is  a  palindromic  polynomial  with  real  coefficients,  i.e.,  the  coefficients  in  front  of  zm~^  and 
zJ  are  real  and  equal  to  each-other  for  j  —  0, . . . ,  m  and  m  >  1: 


..m  |  _  1  |  ~m— 2  i  i  ~  ^^+1  i  i 

2  +  0"m,lZ  +  &m,2^  +  *  *  *  +  CLm^-lZ  2  +  2  + 


+®m,f -l2  2  +  '  '  '  +  Am,22  +  °m,l^  +  1, 


for  m-even, 


I  Ami  = 


(. z  +  l)[zm  1  +  bm!izm  2  H - 1 -bm™-i  ^  1+1+5m  ^iz"1  + 

2  1  +  •  •  •  +  b  m,  2z2  +  bmylz  +  1] ,  for  m-odd. 


(ID 


The  complex  roots  of  a  polynomial  with  real  coefficients  occur  in  conjugate  pairs,  and  if  the 
coefficients  are  palindromic,  then  the  roots  occur  in  inverse  pairs.  Since  the  inverse  pairing  is  not 
necessarily  the  same  as  the  conjugate  pairing,  if  a  root’s  inverse  is  not  its  conjugate,  then  it  must 
be  one  of  a  set  of  four  related  roots  that  satisfy  a  palindromic  quartic  with  real  coefficients.  As  a 
result,  each  palindromic  polynomial  of  even  degree  in  equation  1 1  is  expected  to  factor  into 
quadratic  and  quartic  palindromes  with  real  coefficients.  In  appendix  B,  the  quadratic  and  quartic 
factors  of  |  Am|  for  1  <  m  <  5  are  shown  to  satisfy  the  necessary  and  sufficient  conditions  for  all 
their  roots  to  be  distinct  (complex  conjugate),  and  on  the  unit  circle.  As  a  result,  all  the  zeros  of 
|  Am  |  for  1  <  m  <  5  are  on  the  unit  circle,  and  their  corresponding  angles  are  provided  in 
appendix  D. 

In  this  report,  we  only  consider  m-layered  designs  for  which  all  the  zeros  of  |Am|  are  distinct, 
have  no  multiples,  and  are  on  the  unit  circle.  Based  on  the  results  from  appendix  B,  this  class  of 
designs  is  not  empty.  Under  this  assumption,  Zk  =  eI9k  =  cos  9k  +  X  sin  9k,  and 
z aT1  =  zk  =  e~I6k  =  cos  9k  —  1  sin  9k,  where  zk  is  the  complex  conjugate  of  zk  and  1  <  k  <  |_y J . 
From  the  fact  that 


zk  +  zk  1  =  zk  +  Zk  =  2  COS  6k,  (12) 

a  newly  factored  representation  of  |  Am|  is  obtained,  with  linear  factor  (z  +  1)  for  m-odd,  and  only 
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palindromic  quadratics  with  real  coefficients  involving  the  cosines  of  the  angles  Oi,  92,  ■ . . ,  9 1  ™  , 


Am 


'  LfJ 

[z2  —  2z  cos  9k  +  1]  for  m  even, 

k= 1 

LfJ 

(* + 1)  n  [-~2  —  2z  cos  9k  +  1]  for  m  odd. 

k= 1 


(13) 


Based  on  the  factorization  of  the  determinant  |  Am|  given  in  equation  13,  for  any  m-layered 
design,  if  there  are  m  distinct  roots  on  the  unit  circle,  then  there  are  m  distinct  angles.  The  angles 
{±9k}£=1  correspond  to  the  roots  z  =  e±I9k  for  m  even,  while  the  angles  90  =  n,  {±9k}lkl1 
correspond  to  the  roots  z  =  exe°  =  —  1  and  z  =  c±I(>k ,  respectively,  for  m  odd.  The  [fj 
essential  angles  0  <  9k  <  7T  for  k  =  1.  •  •  • ,  LfJ  were  referred  to  as  the  base  angles  in  Velo  and 
Gazonas  (40).  The  equations  for  the  angles  90  =  7r  and  {9k}k={  in  terms  of  the  design 
parameters  (combinations  of  the  impedance  ratios)  can  be  derived  from  the  condition: 


Dm|  =  0. 


(14) 


The  definition  for  the  matrix  Dm  as  well  as  the  angle-design  parameter  equations  for  up  to  five 
layers  are  included  in  appendices  C-D. 

In  the  rest  of  this  section,  we  prove  that  the  resonance  frequency  spectrum  for  the  discrete  model 

I  —  I 

consists  of  all  the  positive  angles  coterminal  with  0Q  =  tt  or  {±6^}^ ,  as  shown  below 
ulo  =  90  +  210tt  =  (2 10  +  l)7r  (for  m  odd  only), 


<  &h,k  —  @k  +  2Zi7t,  (oi  (i)  h  —  0, 1, 2 . . . ,  (15) 

k  &l2,k  —  ~9k  +  2(h  +  1 ) 7T ,  k  —  1,  .  .  .  ,  LyJ  . 

Equivalently,  the  resonance  frequency  ui  is  expected  to  satisfy  one  of  the  following  relations: 


{coscu  =  —  1  (for  m  odd  only), 
cos  d;  =  cos  9k  for  k  —  1, . . . ,  |_y_|- 


(16) 


Since  the  degree  of  |  Am|  is  m  and  the  degree  of  its  minors  |  Au|  is  m  —  1,  for  i  —  1,  2, ... ,  m, 
the  substitution  of  equation  13  into  equation  9  results  in  the  following  expansion  of  the 
components  xm(i)  of  xm  into  partial  fractions: 
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XmU 


(o  =  z(F{z)-m) 


t?i(-l)i|A1,j| 
|  Am| 


=  z(F(z)  —  f  (0)) 


LfJ 


0  I  V  '■ 
2+1  ~  Zj  z2_2 


o*  fc«+6,* 


ife=l 


2z  cos  0^+ 1 


(17) 


=  (. F(z)-/(0)) 


,  L  2  J 

2+1  2-,  z2_ 

fc=l 


,z2+b. 


i,k 


2z  COS  0fc  +  l 


where  b^ 0  =  0  for  m  even  and  i  —  1, 2, . . . ,  m.  Applying  the  inverse  ^-transform,  one  can  derive 
that 


Z"1 


,  LfJ 

Oj, o~ 

2+1  2~! 


k= 1 


a*  ,  22+6*  ,  2 

2,K  2,fc 

22  — 22  COS  0fc  +  l 


Z’1 


f  tr  g 

2+1 


L  2  J 

+  E 


/  ai5fc-  2(2— cosflfc)  bitk-z  sin 

y  22—  22  COS  0fc  +  l  22— 22  COS  0fc  +  l 


LfJ 

=  6i,0(— l)n  +  a^k  cos  (rc0fe)  +  bi)k  sin  (n0fe), 

fc=i 


(18) 


for  the  following  choices  of  the  coefficients  and  i  —  1,  2, . . . ,  m, 

a*k  cos  0fc  +  b*k 


^2,/c  -  ^2  k  5  (^2  fc  - 


m 


sin  6*; 


for  k  =  1,...,  LyJ. 


(19) 


From  the  fact  that  0  <  0^  <  7r  it  follows  that  sin  6f  7^  0  for  all  k  —  1, ... ,  j+J .  As  a  result,  after 
applying  the  inverse  ^-transform  to  equation  17  and  based  on  equation  18,  a  general 
representation  for  the  stress  terms  in  equation  2  is  obtained, 


Si(n) 


Z  1  (xm(i))  = 


f(n)  * 


LfJ 

bi, o(-l)n  +  E  a,i,k  cos  (■ ndk )  +  bi>k  sin  (n0fc) 

fc=i 


-/(o)  • 


LfJ 

&2,o(-l)n  +  E  cos  Mfc)  +  6i,fc  sin  (■ n0k )  , 

fe=i 


(20) 


where  operation  *  means  convolution,  i  —  1,  2, . . .  m,  and  n  >  0.  The  stress  representation 
previously  derived  in  Velo  and  Gazonas  (40)  can  now  be  recovered  for  the  special  choice  of 
loading  f(n)  =  p  for  n  >  0  (see  appendix  A  for  further  details).  The  stress  representation  in 
equation  20  involves  trigonometric  sums  which  relate  to  the  Chebyshev  polynomials  of  the  first 
and  second  kind  (see  subsection  4.4). 
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In  the  work  that  follows,  we  use  equations  17-18  to  obtain  stress  solutions  when  a  harmonic 
forcing  function  of  the  form  f(n)  =  sin(nu))  with  /( 0)  =  0  and  n  >  0  is  applied,  unless 
mentioned  otherwise.  A  similar  approach  can  be  used  for  a  harmonic  forcing  function  of  the 
form  f(n)  =  cos (nu)  for  n  >  0.  As  discussed  in  Proakis  and  Manolakis  (43),  the  nature  of  time 
(continuous  or  discrete)  is  expected  to  affect  the  nature  of  the  frequency.  The  frequency  for  the 
analog  (continuous)  time  signal  is  related  to  the  frequency  of  the  digital  (discrete)  time  signal 
through  a  linear  transformation  involving  the  time  step  or  the  so-called  sampling  interval.  There 
is  a  linear  transformation  in  frequency 

2  T 

u  =  — oj,  (21) 

m 

from  the  continuous  forcing  function  with  frequency  u  to  the  discrete  forcing  function  with 
frequency  u,  and  time  step 


.  ,  .  .  (  t  2t  \  .  /rL  2r  \  .  (  2 r 

sm(tuj)  =  sm  -5-  •  — uj  — >  sm  •  — uj  =  sin  n  ■  — uj 

V  —  rn  J  \  —  m  )  \  m 

x  m  '  x  m  '  \ 


sin(no;).  (22) 


As  seen  in  figure  1,  t  >  0  represents  the  time  variable,  n  —  \ =  0, 1,  2, . . .  represents  a 

m 

time-related  index,  r  represents  the  wave  transit  time  through  the  layered  strip,  and  m  represents 
the  number  of  layers. 


For  the  choice  of  f(ri)  =  sin(no)),  n  >  0,  unless  mentioned  otherwise,  we  analytically  prove  the 
equivalent  equations  15  and  16,  and  validate  them  through  various  numerical  experiments. 
Substituting  /( 0)  =  0  into  equation  17,  and 

p(z)  =  Z(H  n))  =  Z(s  in(ni,))  =  (,2T2^  +  1)-  <23) 

after  using  equation  18,  it  can  be  shown  that 


Si(n)  =  Z  '(xm(i)) 

2:  sin(d») 


-1 


=  z 


Si(n )  =  bifiZ 


(z2—2z  cosa}+l) 


-1 


z 2  sin(o;) 

{z2— 2z  coscj+1)(2:H-1) 


zsm(u)-(a*kz2+b*kz) 
(z2—2z  cos  9k+l)(z2— 2z  coso)+l) 


or  equivalently, 


Si(n)  =  bi0Z 


-1 


z2  sin(d;) 

(z2—  2z  coso;H-1)(2:H-1) 


if. 

+  E  z 

k=  1 


-1 


zsm(Lj)ja*kz2+b*kz) 
(z2—2z  cos  0k+l)(z2  —  2z  costD+l) 


(24) 


(25) 
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In  the  subsections  that  follow,  we  calculate  the  inverse  ^-transform  of  both  expressions  in 
equation  25  in  order  to  derive  the  stress  formulas  at  resonance  and  non-resonance  frequencies. 

2.1  Stress  Solutions  at  Non-resonance  Frequencies 


Here  we  show  that  the  values  of  the  driving  frequency  d),  which  satisfy  cos  d>  ^  —1  =  cos  6*0  (for 
m  odd)  and  cosd)  /  cos  9  k  for  all  k  =  1, . . . ,  [yj ,  represent  non-resonance  frequency  values.  In 
the  set  of  real  numbers,  this  set  of  frequency  values  is  the  complement  of  the  set  of  values 
identified  in  equation  16.  We  begin  by  finding  the  inverse  c-transform  of  the  two  expressions 
given  in  equations  27  and  29. 


When  cosd;  ^  —1,  the  partial  fraction  expansion 


z 2  sin(d)) 

(z2 -2z  costi+l)(z+l) 


2(l+cosu)) 


2(l+cos  Co) 


implies  that 


7-1 


z 2  sin(u) 


(z2—  2z  coso;+l)(2:+l) 


2(l+coso;) 


z2— 2z  costD+l  z-\- 1 


z(z— cos  a;)  sin  a;  ,  2:sina;(lH-cosa;)  •  ~  z 

z2— 2z  coscIj+1  z2— 2z  cosa;+l 


^rcos(no;)  +  |sin (nu)  —  (— 1) 

uj)  \  /  2  v  /  2(l+cosa;)  v  ' 


When  cos  Cu  ^  cos  9,  the  partial  fraction  expansion, 

z  sin £u(az2  +  bz)  Az(z  —  cosd)) 


+ 


Bzsmu 


(. z 2  —  2 z  cos  u  +  1)  (z2  —  2z  cos  9  +  1)  z2  —  2z  cos  u  +  1  z2  —  2z  cos  u  +  1 

Cz(z  —  cos  9)  D z  sin  9 

z2  —  2z  cos  9  +  1  +  z2  —  2z  cos  9  +  1  ’ 


implies  that 

z-1 


zsmu(az2  +  bz) 


(z2  —  2 z  cos u  +  1  )(z2  —  2zcos9  +  1) 


(26) 


(27) 


(28) 


=  A  cos  (ncu)  +  B  sin  (nuj) 

+C  cos (n9)  +  D  sin(n6)),  (29) 


where  A  —  -C  —  B  =  0(aco~s“+b„„  and  D  =  - Substituting  the 

2(coscj— cos 6)  ’  2(coscj— cos 6/J 7  2smt/(costj— cost/)  ° 

results  from  equations  27  and  29  into  equation  25,  the  following  stress  solutions  for 
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cos  a;  7^  —1  =  cos  90  and  cos  Co  ^  cos  9k  for  all  k  =  1, . . . ,  |_yj  are  obtained: 


s.fn)  =  -  (-1)™  + 

2(l+cos  (i)  1  L)  ^ 


,  LfJ 

°i,0  i  D 

~2~  +  £>i,k 

k=  1 


bifi  sin  a; 
2(l+cosdi) 


LfJ 

+  ^2  A,fc 

fc=i 


cos(ntu) 


LfJ  LfJ 

sin(nw)  +  J2  Q.fccos (n9k)  +  Dijksm(n9k), 


(30) 


fc=i 


fc=i 


where 


a  _  _  ^i,fcsinui  t-j  _  Uj  ^  cos ^  ,fc  cos ;£.)  sinu) 

2(costI>— cosSfc)  ’  2(cosw— cos0fc)  ’  2  sin  6^ (cos  a)— cos  Ok)  '  '  ' 

As  before,  6i>0  =  0  when  the  number  of  layers  m  is  even.  For  a  given  value  of  a),  the  stress 
solutions  in  equation  30,  including  all  the  coefficients  and  angles  9k  for  k  —  1, . . . ,  |_yj ,  can  be 
also  expressed  in  terms  of  the  impedance  ratios  only  (see  equations  6-9,  14,  and  appendices  D 
and  E). 

Given  that  0  <  9k  <  tt  implies  that  sin  9k  ^  0  for  all  k  —  1, ... ,  |_yj  •  This,  together  with  the  fact 
that  cos  a)  f  — 1  and  cos  a;  f  cos  9k  for  all  k  —  1, . . . ,  implies  that  the  coefficients  in  the 
stress  solutions  are  always  bounded  (see  equations  30-31).  Therefore,  for  a  given  m- layered 
design  and  driving  frequency  Co,  the  stress  solutions  in  equations  30  given  by  trigonometric  sums 
of  sines  and  cosines  are  bounded.  However,  as  cos  a;  approaches  —  1  =  cos  9 0  (for  m  odd)  or  any 
of  the  cos  9k  values,  k  =  1, . . . ,  |_™J ,  the  denominators  of  some  coefficients  of  the  stress  solutions 
in  equations  30-31  get  close  to  zero,  allowing  these  solutions  to  take  larger  values  and  approach 
resonance. 

2.2  Stress  Solutions  at  Resonance  Frequencies 

The  stress  solutions  in  equation  30  suggest  that  the  driving  frequencies  Co  identified  in  equation  16 
are  resonance  frequencies.  Here  we  prove  equation  16  and  its  equivalent  form  equation  15,  by 
deriving  the  stress  solutions  at  the  expected  resonance  frequencies,  demonstrating  that  these 
solutions  grow  without  bound  over  time. 

2.2.1  Case  I :  cos  Co  —  —  1  =  cos  00  (for  m  odd  only) 

The  fact  that  cos(cu)  =  —  1  or  equivalently  that  Co  takes  values  u>i0  =  (2 10  +  l)7r  for  l0  —  0, 1,  2  . . ., 
implies  that  f(n)  =  sin(na))  =  sin(n7r)  =  0  for  all  n  >  0.  Asa  result,  due  to  the  zero  input  in 
loading  we  are  unable  to  detect  resonance,  as  confirmed  by  equation  20.  However,  by  choosing 
f(n)  =  cos(no))  =  cos(n7r)  we  overcome  this  obstacle.  Based  on  equations  17-20,  we  expect 
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any  unbounded  stress  terms  to  come  from 


F(z) 


IfJ 


KoZ  a*i,kZ 2  + 

+  2^  y2  _ 


z  +  1 


k= 1 


z 2  —  22  cos  6*^  +  1 


(32) 


after  the  inverse  2-transform  is  applied.  Here  F(z)  =  Z(f(n))  =  Z(cos((mr))  =  which 
after  substituting  into  equation  32,  gives 


Z+1 


LfJ 


Oj,QZ  a*,fc22+^*,fc~ 

2  +  1  _l"  2-,  22  — 2 


2:2— cos  #&  +  ! 


bi,oz2 
(2  +  1)' 


k= 1 
I  —  I 

++  al.22+6*.2 

_l_  »,«  * 

2^  22_2 


bi,oz2 
~  (2  +  1)2 


LfJ 

+  E 

k= 1 


z-(ffl*,fez2+b*,fcz) 

(2:+l)(2:2  —  2z  cos  0&+1) 


k= 1 


2:2— 2z  cos  0^+1 


LfJ 

E 

k= 1 


(33) 


alfe22 +6*^,2 

(2+l)(22  — 22  COS  0fc  +  l)  ' 


Furthermore,  for  the  last  sum  above,  we  have  that 


k= 1 


(2+l)(22 


Z2+ft*,ii+ 

—22:  cos  0fc+l) 


L  2  J  o  ■  /> 

Eui,k  m  _ 2r  sin  _ 

sin^fc  (z+l)(z2—2z  cosO^+l) 

k= 1 


LfJ  * 

Eai,k  #  _ z2  sin  Ok _ 

sin^fc  (2:+l)(2:2— 2z  cos  0^+1) 

k= 1 


LfJ 

_ Z _ 

,k  (z-\-l)(z2—2z  cosOk+l) 

k= 1 


+  Ef-I 


L— J 

I  bi,k  m  z  _J_  -z2  +  (  1+2  cos 0k)z 

2  cos  Ok  z+1  z2  — 2z  cos  0^+1 
k= 1  L 


(34) 


Applying  the  inverse  2-transform  to  the  sums  in  equations  33-34  and  using  equations  18  and  27, 
we  conclude  that  the  only  term  in  equation  33  that  generates  unbounded  stress  values  is  j+jj2  for 
b^o  t -  0,  due  to  the  fact  that 


Z_1 


KoZ2 
(2  +  1)2 


blfi-Wn{-l)  =  b+-l)n{n  +  l). 


(35) 


Here  Wn(y)  represents  the  Chebyshev  polynomial  of  the  second  kind  evaluated  at  y  =  —1.  The 
presence  of  resonance  for  the  choice  of  frequency  values  c ui0  =  (2 l0  +  l)7r,  Z0  =  0,1,2,..., 
proves  that  they  are  part  of  the  resonance  frequency  spectrum  identified  in  equation  15.  Notice 
that  the  resonance  frequency  spectrum  ui0  =  (2 10  +  l)7r  for  l0  =  0, 1,  2, . . .  is  universal  for  all 
the  designs  with  an  odd  number  of  layers  as  it  is  independent  of  any  design  parameters  (e.g., 
impedance  ratios).  This  is  illustrated  in  figure  2  for  the  7-  and  1 1 -layer  case  and  later  on  in 
figure  8  for  the  3-layer  case. 
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(a)  (b) 

Figure  2.  Stress  time  history  for  a  layered  strip  with  impedance  ratios  a\  =  3, 

«2  =  2,  «3  =  1.5,  04  =  2.2,  a 5  =  0.3,  a$  =  1.7,  a-j  =  1.4,  =  3.1, 

«9  =  0.8,  and  aio  =  4  as  applicable;  r  =  1  and  loading 
f(n)  =  cos (nib)  applied  at  £  =  0  for  uj  =  it.  (a)  The  middle  of  the 
second  layer  of  a  7-layered  strip  located  at  £  =  3/14.  (b)  The  middle 
of  the  second  layer  of  an  1 1 -layered  strip  located  at  £  =  3/22. 


2.2.2  Case  II:  COS  6b  —  COS  0fc* 
The  following  expansion 


z  sin  0(az2-\-bz)  _  Az[-2z-\-(l+z2)  cosO -\  .  Bz(z2  —  1)  sin  0  .  CzsinO 

(z2—2z  cos  0-\-l)2  (z2—  2z  cos  0+1)2  J  '  (^2 —2^  cos  0+1) 2  '  22  — 22:  cos  0+1 


(36) 


implies  that 


2  sin  9(az2  +  bz ) 
+2  —  2zcos9  +  l)2 


A  •  ncos(n9)  +  B  ■  nsin(n9)  +  C  sm(n$), 


where  A 


(a  cos  0+6) 

2  sin  0  ’ 


B 


f , and  C 


q+6  cos  0 
2  sin2  6  ' 


2.2.3  Case  Il.a.:  sinco  =  sin0fc* 


(37) 


Applying  equations  27,  29,  and  37  to  equation  25,  the  following  stress  solutions  for  the  frequency 
spectrum  6jtuk*  =  9k *  +  21itt,  l\  =0,1,2,...  are  obtained, 
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Si(n) 


bj,o  sin  dk*  / _ -i  \n  i 

2(l+cos  )  \  / 


6jiOsin0fc* 


LfJ 


0/1  I  /]  T  I  Til 4.1  Z~*  "J-  5  A,‘ 

2(l+cos0^,*)  Z— /  V 

k=l,k^k* 


cos(n9k*) 


LfJ 


sin(n0fc*) 


+  4r  +  nBi^k*  +  C^fc*  + 

k=l,k^k* 

LfJ  LfJ 

+  X]  Ci)fc  cos(ra0fc)  +  A,fcsin(n0fc) 

k=l,k^k*  k=l,k^k* 


(38) 


Here  A^fc,  Ci}k,  A,fc  are  given  from  equation  31  for  a)  =  9k *  +  2/i7t,  /i  =  0, 1,  2, . . while 


A-i  k* 


Kfc*  cosgfc*  +6*|ife0 
2  sin  9k* 


A  L.* 

D  _  l,K  ^  _ 

5  £>i,k*  —  0  5  ^i,k*  — 


+  bU* cos 

2  sin2  9k* 


(39) 


As  before,  bit0  =  0  for  m  even  and  sin  9k  ^  0  for  all  k  —  1, ... ,  |_yj .  The  stress  solutions  in 
equation  38,  including  all  the  coefficients  and  angles  9k  for  A:  =  1, . . . ,  |_^J ,  can  be  also  expressed 
in  terms  of  the  impedance  ratios  only  (see  equations  6-9,  equation  14,  and  appendices  D  and  E). 


All  the  terms  in  the  stress  solution  in  equation  38  are  bounded,  except  for  those  that  are  multiplied 
by  the  time-related  linear  factor  n  given  below 


nAi)k*  cos (n9k*)  +  nBi)k*  sin (n9k*) 


Kfc*  cos 9k*+b*i>k.) 
sin  9k* 

a*k,  cos  [(n  +  1  )9k*] 
sin  9 


cos  (n9k*) 
+  Kk*  cos 


—  a*i  k*  sin (n9k*) 
[( n9k*)]' 


(40) 


where  a*k *  f  0  or  b*k»  f  0  for  i  =  1,2,...,  //?,.  This  implies  that  the  stress  amplitude  is 
expected  to  grow  without  bound  over  time,  proving  resonance  for  the  frequency  spectrum 
=  9k*  +  2Z ! 7r ,  li  =  0, 1,  2, . . .,  included  in  equation  15. 


2.2.4  Case  Il.b.:  sincu  =  —  sin0fc»  =  sin(  —  0k*) 


A  similar  stress  representation  to  equation  38  can  be  obtained  at  the  frequency  values 
uJi2,k*  =  —  9k*  +  2 (/2  +  l)7r,  /2  =  0, 1,  2, . . .,  after  a  sign  adjustment  for  some  of  the  coefficients. 
The  presence  of  resonance  proves  that  this  frequency  spectrum  is  also  part  of  the  resonance 
frequency  spectrum  identified  in  equation  15. 


In  summary,  we  related  to  an  m-layered  Goupillaud-type  strip  with  all  the  zeros  of  |  Am|  on  the 

I  —  I 

unit  circle,  its  corresponding  set  of  angles  90  =  tt  (for  m  odd  only)  and  (±6)fc}^1 .  In  general, 
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the  same  set  of  angles  could  correspond  to  more  than  one  layered  design.  When  the  strip  was 
subjected  to  a  harmonic  forcing  f(n)  =  sin(na))  or  f(n)  =  cos(nio)  at  one  end  and  held  fixed  at 
the  other  end,  we  showed  that  the  positive  angles  coterminal  with  90  =  n  (for  m  odd  only)  and 
{±0k}k=i  represented  values  of  the  driving  frequency  Co,  which  produced  resonance  (see  equation 
15  or  16).  The  complement  of  this  set  of  values  in  the  set  of  real  numbers  was  shown  to  represent 
the  set  of  non-resonance  frequencies. 

In  section  2.3  we  demonstrate  this  approach  and  validate  some  of  the  results  mainly  for  the  two- 
and  three-layer  case.  According  to  appendix  B  and  equations  D-l  and  D-2  in  appendix  D,  a 
Goupillaud-type  strip  with  m  —  2,  3  has  all  the  zeros  of  |  Am|  on  the  unit  circle,  hence  a 
resonance  spectrum  given  by  equation  15.  For  simplicity,  we  continue  to  choose 
f(n)  =  sin(na)),  n  >  0,  unless  mentioned  otherwise. 

2.3  The  Two-layer  Case  (m  —  2) 

For  the  two-layer  case,  it  follows  from  equations  21-22  that  Co  =  tio  and 
f(n)  =  sin (n  ■  tio )  =  sin(na)).  As  a  result,  equation  17  becomes 


zsin(u)) 

2  ( z 2  —  2z  cos  to  +  1) 


where  xi  =  1  +  an  >  1  and  cos^i  =  ^  ±1  from  equation  D-l. 

2.3.1  Stress  Solutions  at  Non-resonance  Frequencies 


(41) 


When  cos  a;  ^  cos  6 1,  we  obtain  the  following  bounded  stress  solutions 


si(n) 

S2(n) 


=  Z  X(x2) 


Xi  (cos  Cj— cos  6\) 


sinwcos(nw)  +  (1  +  cosw)  sin(nw)  —  sinwcos(n6|i)  —  (1+c°^^sln^  siri^G) 


2sinwcos(nw)  +  2  cos  Co  sin(nw)  —  2  sin  Co  cos(n0i)  —  2  c°si^1gS1m  u  sin(n0i) 


(42) 


by  applying  equations  25  and  30  to  equation  41.  As  the  driving  frequency  Co  approaches  any  of 
the  values  from  the  resonance  frequency  spectrum  in  equation  15,  the  denominator 
(cos  a)  —  cos  9i)  of  the  coefficients  of  the  stress  solutions  gets  close  to  zero,  allowing  the  stress 
amplitude  to  take  larger  values. 
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Illustrations  of  the  bounded  stress  solutions  at  non-resonance  frequencies  appear  in  figure  3.  As 
depicted  in  the  plots  in  figure  3  and  discussed  previously  in  this  section,  when  positioned  in  the 
middle  of  the  second  layer,  the  time  interval  at  which  Si(n)  is  reached  is  <  t  <  4ri+^-i 

for  n  >  1  and  i  —  1,2  (see  figure  1).  Beat  phenomena  are  clearly  visible  in  figure  3(b),  as  the 
driving  frequency  approaches  a  resonance  frequency. 


(a)  (b) 

Figure  3.  Stress  time  history  for  a  two-layered  Goupillaud-type  strip  in  the  middle 
of  the  second  layer  located  at  £  =  3/4.  Impedance  ratio  a\  =  1/3, 
angle  6\  =  27r/3,  t  =  1,  and  loading  /(n)  =  sin  (nee)  at  £  =  0.  (a) 

Co  =  7t/4.  (b)  Co  =  0.9  •  6\  =  0.6-7T. 


2.3.2  Stress  Solutions  at  Resonance  Frequencies 

When  cosce  =  cos^i  and  since  =  sin#!,  the  stress  solutions  become 


si(n) 

1 

_  s2(n)  _ 

Xi 

(1fiX1>ncos(,t9i)  + 

2S£ncOS(nel)  + 


n+  jz 


cos  0\ 


sin(n#i 


2  n  +  .  i  „ 

sirr  9 1 


sin(n#i 


(43) 


18 


by  applying  equations  25  and  38  to  equation  41  for  9k »  =  9]_.  Further  simplifications  of  equation 
43  implies 


si(n) 

n 

.  s2{n)  _ 

Xi 

3[(n+l)6 


2  cos[(n+l)#i] 
sin  6\ 


+ 


sin(n6,1) 

Xi 


1  " 

(1  — COS  9\) 

2 

sin2  6>i 


'  cos[(n+|)cos  1(^)j 

VxT 

cosj\n+l)  cos-1  ( ] 

Vxi-i 


r  1  i 

■  n  + 

2 

1 

to 

•  sin 

ncos  - 

Xl 

L  V  xi  /  J 

L  2(X1-1)  J 

(44) 


Given  that  0  <  9i  <  n,  xi  =  (1  +  «i)  and  cos  9i  =  it  follows  that  sin  6\  =  2n/xi  1  and 

xi  Xi 

sin  y  =  ~^=.  Therefore,  all  the  stress  solutions  can  be  expressed  in  terms  of  the  design 
parameter  xi  or  equivalently  in  terms  of  the  impedance  ratio  a ,  only,  as  shown  in  the  last  equality 
of  equation  44. 

As  seen  in  figure  1,  the  stress  values  alternate  between  /(n)  and  Si(n)  in  the  first  layer,  and 
between  si(n )  and  s2(n)  in  the  second  layer,  for  n  =  1,  2, . . ..  The  terms  of  both  stress 
sequences  si(n)  and  s2(n)  become  unbounded  over  time  due  to  the  (time-related)  linear  factor  n 
(see  equation  44).  This  implies  that  the  stress  amplitude  will  grow  without  bound  in  each  layer, 
confirming  resonance.  Based  on  these  results,  we  expect  to  detect  resonance  in  both,  the  first  and 
second  layer,  as  illustrated  in  figure  4.  The  good  agreement  between  the  analytical  results 
generated  from  equation  38,  or  equivalently  equation  44,  and  the  stress  values  generated  from  the 
recursive  equations  2  for  s  i  (n )  and  s2(n),  n  >  1,  is  displayed  in  figure  5.  The  fact  that  the  linear 
factor  n  in  equation  44  multiplies  a  cosine  function  explains  the  sinusoidal  fluctuations  of  the 
stress  solutions  in  figures  4-5. 
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(a)  (b) 

Figure  4.  Stress  time  history  for  a  two-layered  Goupillaud-type  strip  with 

impedance  ratio  a\  =  1/3,  r  =  1,  subjected  to  loading  f(n)  =  sin(nu)) 
with  Co  =  9\  =  27r/3.  (a)  Middle  of  the  first  layer  located  at  £  =  1/4. 
(b)  Middle  of  the  second  layer  located  at  £  =  3/4. 


Figure  5.  The  values  of  the  two  stress  sequences  (a)  si(n),  (b)  S2(n),  obtained 
from  the  recursive  equations  2  and  the  analytical  equation  38,  for  the 
two-layered  Goupillaud-type  strip  with  impedance  ratio  a\  =  1/3, 
r  =  1,  subjected  to  loading  /(n)  =  sin(nd;)  with  Co  =  9\  =  27t/3. 
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2.4  The  Three-layer  Case  (m  =  3) 


For  the  three-layer  case,  it  follows  from  equations  21-22  that  and 

f(n)  =  sin (n  ■  =  sin (nu).  As  a  result,  the  first  equation  of  17  becomes 

2(1+02)  z2 +4(1— 0.2)2+2(1+02) 

X2(z+1)(z2-2zcos91+1) 

_ J+ _  (A 

X2  (z2 — 2z  cos  0i+l)  ’  ' 

_ 8 z^_ _ 

X2(z+1)(z2-2z  cos  0i+1) 

where  X2  —  (1  +  «i)(l  +  ct2)  >  1,  cos  90  =  —1,  and  cos$i  =  X^L  7^  ±1  from  equation  D-2  in 
appendix  D.  Applying  in  equation  45  partial  fraction  expansion  as  shown  and  equation  25, 
results  in 


z  sin (u) 


x3  = 


( z 2  —  2  z  cos  Cj  +  1) 


si(n) 

2 

2ol2 

,  [(1+02)  cos  01  +  (1— a2)](z+l) 

" 

1+cos  6 1 

z+1 

z2  —  2z  cos  0i  +  l 

S2i.11) 

-  Z~1(xa)  -  1  •  22sin(<i) 

A.z 

Z  lX3j  z  lX2  (22_2zcosii+l) 

z2—2z  cos  0i  +  1 

_  s3{n ) 

4 

1  .  (1+2  cos  Q\)z— 1 

1+COS  61 

2+1  '  z2  —  22  cos  9 i+l 

(46) 
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2.4.1  Stress  Solutions  at  Non- resonance  Frequencies 


When  cos  a;  ^  —1  =  cos  90  and  cosd)  7^  cos  9 1,  bounded  stress  solutions  are  obtained  by  applying 
equation  30  to  equation  46  and  further  simplifying  the  results, 


si(n) 

—  2c»!2  sino; 

[(l+a2)  cosa;+(l— q2)]  sin  a; 

( 1 +cos  uj  )  ( 1 +cos  9 1 ) 

(cos  a;— cos  9\ )  (1+cos  to) 

s2+) 

=  — { 

X2  L 

0 

(-1)”+ 

2  sin  a) 

(cos  lj— cos  6 1) 

cos  (nd)) 

.  s3(n)  _ 

2  sin  Cj 

2(2  cos  cD+l)  sin  a) 

(l+cosd;)(l+cos  6 1) 

(cos  a;— cos  9\ )  (1+cos  uj) 

(I+CK2)  coscj+(1— a2) 
(cos  uj— cos  Qi  ) 


[(l+q2)  cos  0i+(l— 0^2)]  sind) 
(cos  uj— cos  0i)  (1+cos  6\) 


+ 


2  cos  uj 

(cost!;— cos  0i) 


sin  (/id;)  — 


2  sind; 

(cos  uj— cos  0\ ) 


cos(n6+  (47) 


2(2  cosd;— 1) 
(cos  Cj— cos  6\) 


2(2  cos 0i+l)  sino; 
(cos  a;— cos  0i)  (1+cos  6\) 


[(l+a2)  cos  0j+(l— q2)]  sin  a; 
(cos  a)— cos  0 1 )  sin  0\ 


2  cos  0 1  sin  Cj 
(cos  a;— cos  6\ )  sin  0\ 


sin(n0i)  }. 


2(2  cos  0i  —  l)  sin  a; 
(cos  a;— cos  0i )  sin  0i 


As  the  driving  frequency  Cj  approaches  any  of  the  values  from  the  resonance  frequency  spectrum 
in  equation  15,  the  denominators  (cosd)  +  1)  or  (cosd)  —  cos  6 \ )  of  the  coefficients  of  the  stress 
solutions  in  equation  47  get  close  to  zero,  allowing  the  stress  amplitude  to  take  larger  values. 


Illustrations  of  the  bounded  stress  solutions  at  non-resonance  frequencies  are  given  in  figure  6. 

As  depicted  in  figure  6  and  discussed  previously  in  this  section,  when  positioned  in  the  middle  of 
the  second  or  third  layer,  the  time  interval  at  which  s*(n)  is  reached  is  4n+^~3  <  t  <  4n+^~1,  as 
applicable,  for  n  >  1  and  i  =  1,  2,  3  (see  figure  1).  The  good  agreement  between  the  analytical 
results  generated  from  equation  30,  or  equivalently  equation  47,  and  the  stress  values  generated 
from  the  recursive  equations  2  for  si(n),  n  >  1,  as  well  as  the  error  introduced  by  the 
computational  approximation  of  quantities  (such  as  sine  and  cosine  values),  are  displayed  in 
figure  7. 
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(a)  (b) 

Figure  6.  Stress  time  history  for  a  three-layered  Goupillaud-type  strip  with 

=  3,  a2  =  2(1— cos(V4))  “  !»  #i  =  f »  t  =  1.  subjected  to  loading 
/(n)  =  sin(nu))  with  <Z>  =  (a)  Middle  of  the  second  layer  located  at 

£  =  1/2.  (b)  Middle  of  the  third  layer  located  at  £  =  5/6. 


Figure  7.  The  three-layered  Goupillaud-type  strip  with  an  =  3, 

“2  =  2(1— cos(7t/4))  -  b  6>i  =  f ,  r  =  1,  subjected  to  loading 
/(n)  =  sin(nu))  with  u  =  (a)  The  values  of  the  stress  sequence 

si(n),  obtained  from  the  recursive  equations  2  and  the  analytical 
equation  30.  (b)  The  difference  of  the  expressions  for  si(n)  obtained 
from  the  recursive  equations  2  and  the  analytical  equation  30. 
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2.4.2  Stress  Solutions  at  Resonance  Frequencies 


2.4.2.1  Case  I.  cosw  =  —  1  =  cos0o 

The  stress  formulas  in  the  ;j- space  for  the  three-layer  case  in  equation  45,  indicate  that  the 
sequences  Si(n)  and  s3(n)  increase  without  bound  because  cos(cD)  =  — 1  =  cos  90  corresponds  to 
a  zero  of  the  denominator  z  =  exe°  —  — 1.  Alternatively,  we  do  not  expect  unbounded  growth  for 
the  sequence  s2(n),  since  z  =  exe°  =  —  1  is  not  a  zero  of  the  denominator.  Figure  8  illustrates 
this  behavior  for  each  of  the  stress  sequences.  Similar  behavior  was  observed  in  the  numerical 
experiments  for  the  seven-  and  eleven-layer  case  (see  illustration  in  figure  2)  where  the  stress 
sequences  with  odd  subindices  increase  without  bound  while  the  stress  sequences  with  even 
subindices  are  bounded. 


(a)  (b) 

Figure  8.  Stress  time  history  for  a  three-layered  Goupillaud-type  strip  with 
«i  =  3,  a2  =  2(1—008^/4))  ~  T  =  !>  subjected  to  loading 
f(n)  =  cos (nui)  with  Co  =  n.  (a)  Middle  of  the  second  layer  located  at 
£  =  1/2.  (b)  Middle  of  the  third  layer  located  at  £  =  5/6. 


2.4.2.2  Case  Il.a. 

When  cos Cj  =  cos  9 1  and  sin ui  =  sin 9\,  and  by  applying  equation  38  to  equation  46  the  stress 
solutions  become, 
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si(n) 


Oi  2 


ex. 2  sin  0i  _  [(1+0:2)  cos  0i  +  (l— <32)] 

(1+cos  0i)2  2  sin  0i  1 


S2{n) 


2  sin  0i 

X2(l+COS0i)2 


0 


(-1  r 


2 

X2 


sin  t/i 


cos(n0i) 


S3W 


sin  Oi _ (2  cos  0i~  1) 

(1+cos  0i)2  sin  0i 


"  (1— 02)  cos  6>i  +  (1+q2)  ,  [(1+o2)  cos  0j  +  (l  —  o2)]  " 

2  sin2  0i  '  2(l+cos0i)  !> 


+ 


2 

X2 


n  + 


1 

sin2  0i 


2— cos  0i  _i_  (2  cos  0i  +  l) 
sin2  0i  1+cos  0i 


sin(n0i). 


(48) 


Since  the  focus  here  is  the  study  of  resonance,  the  unbounded  part  of  each  of  the  stress  solutions 
containing  the  time-related  factor  n  can  be  simplified  as  follows: 


[(l+«2)  cosfli  +  (l-o;2)] 

2  sin  0i 


COS+$i)  + 


[(I+02)  COS01+(1  — 02)] 
2(l+cos0i) 


sin(n0i) 


2  n 
X2 


^COS^)  +801+00 


(2  cos  01  —  1) 
sin  0i 


COS+0i)  + 


(2  cos  0i+l) 
1+COS  01 


sin +00 


[(l+o2)  cos  01  +  (1— o2)] 
2  cos  sin  0i 


cos  [(n  +  |)  Oi] 


20ly/X2 

(l+Of)(x2-l) 


COS 


X2-2 

X2 


2 n 
X2 


cos[(n+l)0i] 
sin  0i 


cos[(n+l)cos  1  ++] 


cos[(n,+  §)  + 

0  -\  •  n 

cos  ~y  sin  0i 


VX2 

(X2-1) 


•  cos 


*2-2 

X2 


(49) 


Given  that  0  <  91  <  n,  \2  =  (1  +  Q?i)(l  +  a2),  and  cos0i  =  it  follows  that 

sin  0i  =  2^~1-  Therefore,  the  stress  solutions  can  be  expressed  in  terms  of  the  impedance 

ratios  an  and  a2  only.  This  is  demonstrated  in  the  last  equality  (equation  49). 
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As  seen  in  figure  1,  the  stress  values  alternate  between  f(n)  and  si(n)  in  the  first  layer,  between 
si(n)  and  .sv (n)  in  the  second  layer,  and  between  s^in)  and  s3(n)  in  the  third  layer,  for 
n  —  1,2,....  All  three  stress  sequences/solutions  si(n),  S2(n),  and  33(71)  become  unbounded 
over  time  due  to  the  (time-related)  linear  factor  n.  This  implies  that  the  stress  amplitude  will 
grow  without  bound  in  all  three  layers  (figure  9).  The  good  agreement  between  the  analytical 
results  generated  from  equation  38,  or  equivalently  equation  48,  and  the  stress  values  generated 
from  the  recursive  equations  2  for  si(n),  n  >  1,  as  well  as  the  error  introduced  by  the 
computational  approximation  of  quantities  (such  as  sine  and  cosine  values)  are  displayed  in  figure 
10.  The  fact  that  the  linear  factor  n  multiplies  a  cosine  function  as  shown  in  equation  49, 
explains  the  sinusoidal  fluctuations  of  the  stress  solutions  in  figures  9-10. 


(a)  (b) 

Figure  9.  Stress  time  history  for  a  three-layered  Goupillaud-type  strip  with 
=  3,  a2  =  2(i-cos(7r/4))  ~Ur  =  1,  subjected  to  loading 
f(n)  =  sin(no))  with  u  =  9\  =  (a)  Middle  of  the  second  layer 

located  at  £  =  1/2.  (b)  Middle  of  the  third  layer  located  at  /  =  5/6. 
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Figure  10.  (a)  The  values  of  the  stress  sequence  si(n),  obtained  from  the 
recursive  equations  2  and  the  analytical  equation  38.  (b)  The 
difference  of  the  expressions  for  si(n)  obtained  from  the  recursive 
equations  2  and  the  analytical  equation  38  for  a  three-layered 
Goupillaud-type  strip  with  on  =  3 ,  a2  =  2(i-cos(7r/4))  ~ *  1 2>  r  =  1- 
subjected  to  loading  /(n)  =  sin(nu;)  with  uj  =  8\  =  j. 


3.  Analytical  Formulas  for  the  Natural  Frequency  Spectrum  of  a 
Free-fixed  Goupillaud-type  Layered  Elastic  Strip 


In  this  section,  we  derive  the  natural  frequencies  of  a  free-fixed  m-layered  Goupillaud-type  strip 
in  two  different  ways: 

1.  By  simplifying  and  then  solving  the  frequency  equation  for  2  <  m  <  5. 

2.  By  extending  the  results  of  the  resonance  frequency  spectrum  from  the  discrete  model  to 
the  continuous  model  for  any  m  >  2.  In  the  continuous  model,  the  strip  is  subjected  to  a 
continuous  harmonic  forcing  function  at  one  end  and  held  fixed  at  the  other  end. 

3.1  The  Frequency  Equation:  Explicit  Formula  for  the  Natural  Frequency  Spectrum 

The  frequency  equation  for  a  free-fixed  m-layered  Goupillaud-type  strip  is  derived  in  this  section, 
by  generalizing  the  procedure  described  in  Graff  (33).  The  frequency  equation  is  simplified 
through  a  transformation  of  the  spatial  variable,  to  derive  the  natural  frequency  spectrum  for 

2  <  m  <  5. 
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The  governing  equations  for  the  displacement  u(x,  t )  in  each  layer  of  the  strip  are: 


d2Uj(x,t)  _  2  d2Uj(x,t) 

dt2  ~  Li  dx2 


for  Xi_ i  <  x  <  Xi, 


i  —  1, . . . ,  m, 


(50) 


where  0  =  x0  <  Xi  <  ■  ■  ■  Xi  <  ■  ■  ■  <  xm  =  L,  L  represents  the  length  of  the  strip,  c,  represents 
the  wave  speed  in  the  zth  layer,  x  represents  the  spatial  variable,  and  1  >  0  represents  the  time 
variable.  Seeking  harmonic  solutions,  we  may  represent  u(x,  t )  =  U(x)u>(t)  and  the  spatial 
component(s)  of  the  displacement  for  each  layer  as 


Ui(x)  —  Ai  smj3iX  +  Bi  cos/3jX  for  Xj_i  <  x  <  x^  i  —  (51) 


where  /%  =  cu/c*  for  i  =  1 , ,m,  and  u>  represents  the  frequency  of  the  harmonic  motion. 
Similarly,  -0(f)  can  be  represented  in  the  form  -0(f)  =  A*  sin  cut  +  B*  cos  cut  or  in  exponential 
form  as  0(f)  =  A**eXujt  +  B**e~XuJt,  where  X  =  \/—\.  As  a  result,  after  applying  the  free-fixed 
boundary  conditions  at  x  =  0  and  x  =  L,  respectively, 


dUi(0) 

dx 


Um(L )  —  0, 


(52) 


and  the  continuity  of  stress  and  displacement  at  each  layer  interface,  results  in  the  following  linear 
system  of  equations: 


Ai  =  0, 


for  i  =  1, . . . ,  m  —  1, 


Ai  sin  (f,  x^j  +  Bi  cos  =  Ai+1  sin  +  Bi+1  cos 

oti  Ai  cos  -  Bi  sin  =  Ai+1  cos  -  Bi+i  sin  ( -^xij  for  i  =  1, . . . ,  m  -  1, 


(53) 


Am.  sin  A- 


+  Bm  cos  (^TXm)  = 


The  relations  a*  =  -  +c^+1  =  f1  •  were  applied  equation  53,  using  the  fact  that  c  -  p  =  E/c. 
The  condition  that  the  determinant  of  the  system  matrix  for  equation  53  must  vanish  generates  the 
frequency  equation  for  to.  Solving  the  frequency  equation  and  obtaining  a  general  formula  could 
be  a  challenging  task.  However,  if  the  spatial  variable  x  of  the  wave  equation  is  replaced  with  the 
new  variable  £  =  J(Jl  ^0.  the  frequency  equation  can  be  simplified.  This  transformation  converts 
the  problem  to  a  Goupillaud-type  layered  strip  with  the  same  length  L  as  the  wave  transit  time  r 
through  the  strip  (X  =  r)  (see  Velo  and  Gazonas  (40)  and  section  2).  Under  this  transformation, 
the  layers  have  equal  lengths  of  0  and  wave  speed  of  unity.  In  addition,  the  previous  coordinates 
x0  =  0,  Xj,  and  xrn  =  L  arc  replaced  by  the  new  coordinates  £0  =  0,  0  =  0,  and  =  r  for 
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j  =  1,2 , ,m,  while  the  impedance  ratio  a*  =  =  X EiPi  --  remains  the  same  for 

^ 'Ji  + 1  Y  -Ei  +  lPi  +  l 

i  =  1,  2, . . . ,  m  —  1.  As  before,  E,,  p{  represent  the  actual  material  properties  for  the  physical 
case,  while  Ei,  pi  represent  the  material  properties  of  the  / 1 h  layer  for  the  simplified  case,  with 
Ei  =  pi  =  yjEipi .  As  seen  next,  both  designs  of  the  strip  share  the  same  natural  frequency 
spectrum  under  free-fixed  boundary  conditions. 

As  a  result,  equation  53  can  be  transformed  in  a  new  matrix-vector  form  as  follows: 


Amv  =  6, 


(54) 


where 


A 


m 


1  0  0  0 
sin 


ot\  cos  71 

0 

0 
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,v  = 
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(55) 
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—  sin  7i 

—  cos  7i 
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—  cos  7i 
sin7i 
cos  72 
-OL2  sin  72 


0 

0 

—  sin  72 

—  cos  72 
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2m  x  1 


2mx2m 

0  0  0 

—  COS  72 

sin  72 


0 

0 

0 
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0 

0 
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sin7m_i 
OLm  —  l  COS  7m  —  1 
0 


COS  7m -1 
1  sin  7m  _1 
0 


0 


-  sin  7m- 1 

-  COS  7m- 1 


sin  7m 


0 

-  COS  7m- 1 
sin  7m -i 
COS  7m 


(56) 
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Here  7 \  =  ^o;  =  —  for  i  —  1, . . . ,  m.  Seeking  a  nonzero  solution  for  equation  54,  the 
determinant  of  the  system  matrix  Am  must  vanish: 


Am|  =  0. 


(57) 


Equation  57  gives  the  transformed/simplified  frequency  equation  for  u>.  Equation  57  can  be 
successfully  solved,  using  symbolic  algebra  software,  for  a  strip  with  up  to  five  layers  (see 
appendix  F).  Solving  the  frequency  equation  57  for  the  general  m-layer  case,  poses  a 
computational  challenge. 

Based  on  these  derivations,  a  comparison  of  equation  14  with  equation  57,  equation  D-l  with 
equation  F-l,  equation  D-2  with  equation  F-3,  equation  D-4  with  equation  F-5,  and  equation  D-6 
with  equation  F-9  reveals  that  the  two  matrices  Dm  and  Am,  defined  in  equation  C-l  and 
equations  55-56,  have  the  same  dimensions  and  cosine  roots  of  their  determinants  subject  to  the 
transformation 


cos  8  =  cos  — u. 

m 


(58) 


Also  notice  that  the  pairs  of  the  frequency  equations  F-l,  F-3  and  F-5,  F-9,  have  a  similar  pattern 
to  equations  D-l,  D-2  and  D-4,  D-6,  respectively.  Equations  58  can  be  justified  from  equations 
21-22  and  the  fact  that  the  natural  frequencies  of  a  free-fixed  elastic  strip  correspond  to  the 
resonance  frequencies  when  a  (continuous)  forcing,  which  varies  harmonically  with  time  is 
applied  at  the  free  end  of  the  strip. 

3.2  Heuristic  Approach:  Implicit  General  Formula  for  the  Natural  Frequency  Spectrum 

Using  the  linear  frequency  transformation  in  equation  21,  the  resonance  frequency  spectrum  for 
the  continuous  model  (where  the  strip  is  subjected  to  a  continuous  harmonic  forcing  at  one  end 
and  held  fixed  at  the  other  end)  can  be  recovered  from  the  resonance  frequency  results  in  equation 
15  for  the  discrete  model.  This  results  in  the  following  implicit  formula  for  the  natural  frequency 
spectrum  of  a  free-fixed  m-layered  Goupillaud-type  strip: 

■  [Bo  +  2Z07r]  =  .  (2 10  +  1)7 r,  (for  m  odd  only), 


<^h,k  —  ■  [Bk  +  2Zi 7r]  —  ’  [8k  +  2Zi7t], 


lo,  h,l 2  —  0, 1, 2, ... , 


(59) 
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Substitution  of  the  angle  values  from  equations  D-l,  D-3,  D-5,  and  D-7  into  equation  59,  results 
in  the  natural  frequency  spectrum  for  the  free-fixed  Goupillaud-type  layered  strip  with  up  to  five 
layers,  equations  F-2,  F-4,  and  F-6,  derived  independently  in  appendix  F  from  the  frequency 
equation  57.  Although  a  rigorous  proof  for  the  equation  59  is  not  provided,  our  theoretical  and 
numerical  results  validate  it  as  the  general  natural  frequency  spectrum  of  a  free-fixed  m-layered 
Goupillaud-type  strip  with  the  zeros  of  the  corresponding  matrix  |Am|  on  the  unit  circle  (see 
matrix  Am  defined  in  the  z-space  in  equation  7).  The  zeros  of  |  Am|  on  the  unit  circle  generate 
the  respective  angles  90  =  tt  (for  m  odd  only)  and  needed  to  recover  the  natural 

frequency  spectrum. 

3.3  Interpretation  and  Illustration  of  the  Natural  Frequency  Results 

The  explicit  natural  frequency  equations  (F-2,  F-4,  and  F-6)  derived  from  the  frequency  equation 
57,  as  well  as  the  general  frequency  spectrum  derived  from  the  implicit  equation  59,  suggest  that 
for  a  free-fixed  m-layered  Goupillaud-type  strip,  the  natural  frequency  spectrum  does  not  depend 
on  the  actual  length  of  the  strip.  According  to  these  formulas,  the  natural  frequency  spectrum 
depends  only  on  certain  combinations  of  the  impedance  ratios  an,  a2, . . . ,  am- 1  and  (it  is 
inversely  proportional  to)  the  (equal)  wave  travel  time  —  for  each  layer. 

The  following  illustration  validates  these  expectations  and  demonstrates  the  simplicity  of 
predicting  the  natural  frequency  spectrum  when  using  the  general  equation  59,  or  solving  the 
simplified  frequency  equation  57,  versus  solving  directly  the  frequency  equation  generated  from 
equation  53.  Our  numerical  experiments  with  up  to  five  layers  have  shown  that  the  natural 
frequency  spectrum  can  be  predicted  using  either  equation  59  and  appendix  D,  or  appendix  F, 
even  when  a  symbolic  algebra  software  such  as  Maple  cannot  solve  the  frequency  equation 
obtained  from  equation  53. 

3.3.1  Case  a.  Two-layered  Goupillaud-type  Elastic  Strip  with  Unequal  Layer  Lengths 

•  Layer  lengths:  L1  =  \/2,  L2  =  3\/2 

•  Length  of  the  strip:  L  =  Li  +  L2  —  4  a/2 

•  Wave  speed  per  layer:  c\  —  1,  c2  =  3 

•  Elastic  modulus  per  layer:  E1  —  2  +  \/3,  E2  =  3(2  —  v/3) 

•  Wave  travel  time  per  layer:  —  =  —  =  —  =  \[2 

r  J  ci  C2  m  v 

•  Transit  time  through  the  strip:  r  =  2y/2 

•  Impedance  ratio:  an  =  ^ 
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•  Frequency  equation  from  equation  53: 


sin  (y/2u) 
Q?i  cos  (a/2  cu) 
0 


cos 


0 

(V2u) 

-ai  sin  {y/2  cu) 
0 


sill  I  ^  cu 


—  cos  ^  cu 
sin  (^0-  cu 


cos  |  4^  cu 
sin  ^  cu 
cos  ^  cu 


=  0. 


(60) 


3.3.2  Case  b.  Simplified  Two-layered  Goupillaud-type  Elastic  Strip  with  Equal  Layer 
Lengths 

This  free-fixed  strip  has  the  same  natural  frequency  spectrum  as  the  one  in  Case  a. 

•  Layer  lengths:  Lx  =  =  L2  =  §  =  \  =  \/2 

•  Length  of  the  strip:  L  =  Lx  +  L2  =  r  =  2y/2 

•  Wave  speed  per  layer:  c\  —  c2  =  1 

•  Elastic  modulus  per  layer:  E\  =  2  +  y/3,  E2  =  2  —  \/3 

•  Wave  travel  time  per  layer:  4^  =  4^-  =  |  =  \/2 

•  Transit  time  through  the  strip:  r  =  2\/2 

•  Impedance  ratio:  an  =  |/  =  =  an 

•  Frequency  equation  from  equation  57: 


sin  (y/2u) 
a\  cos  (v/2  cu) 
0 


0  0  0 
cos  (\/2w)  —  sin  (a/2  cu)  —  cos  (y/2u) 

— a\  sin  (\/2  cu)  —  cos(\/2cu)  sin  (v^  cu) 

0  sin  (2  y/2  cu)  cos  (2  a/2  cu) 


=  0. 


(61) 


3.3.3  The  Predicted  Natural  Frequency  Spectrum 

The  common  natural  frequency  spectrum  for  the  free-fixed  strip,  for  Cases  a.  and  b. 


“  2^2  (e  +  2^l7F)  “  l//2  +  n/2^1’ 


12V2  V2 


h,l2  —  0, 1,  2, ... , 


“  573  (_l  +  2^2  +  1)7r)  “  S 


(62) 
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can  be  derived  easily  and  accurately  from  equation  59  and  equation  D-l  for 
9i  =  cos-1  )  =  §,  °r  equation  F-2.  The  Maple  software  recovers  most  of  these  results  by 
finding  the  zeros  of  the  determinants  previously  given  in  Cases  a.  and  b.  As  expected,  the 
frequency  equation  in  the  (simplified)  Case  b.  is  much  easier  to  solve  in  Maple  than  the  original 
frequency  equation  in  Case  a. 

3.3.4  Case  c.  Design  Modification  That  Gives  a  Desired  Frequency  Spectrum  Within 
Limitations 

One  way  to  modify  the  natural  frequency  spectrum  equation  62  is  to  modify  any  of  the  two 
layered  strips  given  in  Cases  a.  and  b.,  with  xi  —  (1  +  «i)  and  ()\  =  |,  by  adding  a  new  layer  in 
front  of  the  existing  two  layers.  The  length  of  the  new  layer  must  be  chosen  so  that  the  strip 
remains  of  Goupillaud-type.  The  new  three-layered  strip  is  characterized  by 
X2  —  O-  +  «i)(l  +  ol 2)  with  a\  to  be  determined  and  a*2  =  ot\.  For  a  desirable  value  of  the 
cosine  of  the  base  angle  0*  given  by  V  =  cos  0*  =  one  can  derive  the  necessary  formula  for 
the  unknown  impedance  ratio  a{: 


cos  6*i  <  V  =  cos  9\  <  1, 


(63) 


using  the  fact  that  y2  =  (1  +  ®i)Xi-  For  instance,  for  the  choice  of  0\  =  |  we  have  that 
cos  9 1  =  cos  (|)  <  V  =  cos  9\  =  cos  (|) ,  which  results  in  the  natural  frequency  spectrum: 


( 


(64) 


for  the  free-fixed  three-layered  strip  with  impedance  ratios  a\  =  —1  + 
and  the  (equal)  wave  travel  time  per  layer  ^  |  —  \/2. 


2(  1— cos 


2-V 
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4.  Other  Applications  of  the  Frequency  Results 


The  properties  of  the  materials  used  in  this  section  appear  in  appendix  G.  Using  the  definition  of 
the  wave  speed  c  =  y/E/ p,  the  characteristic  impedance  pc  in  each  layer  is  calculated  in  terms  of 
the  material  properties  as  y/Ep.  Once  the  length  of  one  layer  is  chosen,  the  length  of  each 
remaining  layer  is  determined  according  to  the  material  properties  to  provide  equal  wave  travel 
time:  E  =  T-^  or  Li+1  =  \J  •  ff  L{.  Here  Lv  cv  E:j  and  pn  for  j  —  i,  i  +  1,  represent  the 
length,  wave  speed,  elastic  modulus,  and  material  density  for  the  Ah  and  (i  +  l)th  layer, 
respectively,  i  =  1, ...  pm  —  1. 

4.1  One-dimensional  Layered  Media  with  a  Common  Frequency  Spectrum 

According  to  equations  14  and  15,  the  resonance  frequency  spectrum  for  the  discrete  model 
depends  only  on  certain  combinations  of  the  impedance  ratios  cci,  a2,  •  •  • ,  am_ 1-  Such 
combinations  of  impedance  ratios  when  2  <  m  <  5  are  represented  by  the  design  parameters 
Xm-i  for  m  =  2,  3,  and  Xm-i>  Tm-i  for  m  =  4,  5.  There  are  (theoretically)  infinitely  many 
m-layered  designs  of  a  Goupillaud-type  elastic  strip,  which  have  the  same  values  of  the  design 
parameters  Xm- i,  Tm-i  for  2  <  m  <  5,  hence  the  same  resonance  frequency  spectrum  for  the 
discrete  model. 

As  for  the  natural  frequency  spectrum  of  a  free-fixed  m-layered  Goupillaud-type  strip  based  on 
equation  59,  it  depends  not  only  on  combinations  of  the  impedance  ratios  cti,  a2,  ■  ■  ■ ,  ctm- i  but 
also  on  the  (equal)  wave  travel  time  ^  for  each  layer.  Therefore,  to  obtain  designs  with  a 
common  natural  frequency  spectrum,  besides  having  the  same  values  of  the  relevant  design 
parameters,  the  layer  lengths  have  to  be  chosen  so  that  they  all  have  the  same  wave  travel  time. 
Another  way  to  find  designs  with  a  common  natural  frequency  spectrum  was  previously 
illustrated  in  subsection  3.3  using  a  transformation  of  the  spatial  variable. 

For  the  three-layer  case  where  Xi  =  (1  +  o , )  ( 1  +  a2),  the  tungsten-lead-aluminum  configuration 
with  «  5.828  and  a2  ~  0.952  has  the  same  value  of  X2  ~  13.33  as  the  lead-aluminum-copper 
configuration  with  a.\  ~  0.952  and  a2  ~  5.828.  Another  material  design  with  y2  ~  13.33  is  the 
aluminum-steel-brass  configuration  with  ~  0.331  and  a2  ~  9.016.  In  addition,  in  order  for 
these  material  designs  to  have  a  common  natural  frequency  spectrum  for  the  free-fixed  strip,  only 
the  length  of  one  layer  in  one  of  the  designs  may  be  chosen  at  random.  The  lengths  of  the 
remaining  layers  in  all  the  three  designs  are  then  determined  to  provide  equal  wave  travel  time. 
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For  the  four-layer  case,  we  have  that  %3  =  (1  +  oii)(l  +  a2)(l  +  a3)  and  V3  =  a±a3  —  1.  The 
aluminum-(glass  fiber  reinforced  material)-nickel-bronze  configuration  with  an  ~  2.27, 
a2  ~  0.698  and  a3  —  1/an  tv  .44,  has  almost  the  same  value  of  %3  =  8  and  T3  =  0,  as  the 
homogeneous  strip  (all  “four  layers”  are  made  of  the  same  elastic  material).  As  a  result,  these 
designs  share  the  same  resonance  frequency  spectrum  when  subjected  to  a  discrete  harmonic 
forcing  at  one  end  and  held  fixed  at  the  other  end.  The  common  frequency  spectrum  can  be 
derived  from  equation  15  after  substituting  y3  =  8  and  T3  =  0  into  equation  D-5.  In  addition,  a 
common  wave  travel  time  for  each  of  the  four  layers  in  both  designs  would  ensure  the  same 
natural  frequency  spectrum  under  free-fixed  boundary  conditions. 

Furthermore,  according  to  equation  D-5  and  table  E-l,  these  designs  with  parameters  X3  =  8  and 
r3  =  0,  and  base  angle  values  91  <opt  =  7t/4  and  92)0pt  =  37t/4  obtained  for  j  =  2,  are  expected  to 
be  optimal,  i.e.,  provide  a  stress  amplitude  of  double  the  loading  when  a  Heaviside  step  in  stress 
loading  is  applied  at  one  end  of  the  strip  while  the  other  end  is  held  fixed. 

4.2  Resonance  Frequencies  and  Optimal  Designs 

For  the  m-layered  strip  subjected  to  a  Heaviside  loading  (40),  it  was  shown  that  when  the  optimal 
values  of  the  base  angles  9k:0pt  for  k  —  1, ... ,  [yj ,  given  in  appendix  E,  were  substituted  into  the 
respective  angle-design  parameter  equations  in  appendix  D,  optimal  m-layered  designs  were 
generated  for  2  <  m  <  5.  These  optimal  designs  provide  the  smallest  stress  amplitude  of  double 
the  loading  in  all  layers,  for  all  times. 

However,  when  these  optimal  angle  values  are  substituted  into  equations  15  and  59,  they  generate 
the  resonance  frequency  spectrum  for  the  discrete  and  continuous  model,  respectively.  This 
means  that  for  the  optimal  designs  identified  in  appendices  D-E,  the  displacement  and  stress 
amplitude  will  grow  without  bound  over  time  when  a  harmonic  forcing  with  a  resonance 
frequency  is  applied  instead  of  the  Heaviside  loading.  As  a  result,  for  these  designs,  depending 
on  the  type  of  loading,  one  can  get  either  the  best  or  the  worst  outcome  when  it  comes  to 
controlling  the  stress  amplitude  (see  the  illustration  in  figure  1 1).  As  previously  seen  in  figure  1, 
the  stress  values  in  a  two-layered  strip  are  given  by  the  sequences  f(n),  si(n),  and  s2(n),  n  >  0. 
For  the  example  considered  in  figure  1 1(a),  we  have  that  |  f(n)  \  <  1  for  all  n  >  0,  therefore  the 
stress  amplitude,  which  exceeds  unity,  is  determined  by  the  other  two  stress  sequences. 

Similarly,  for  the  example  considered  in  figure  1 1(b)  we  have  that  f(n)  =  1  for  all  n  >  0, 
therefore  the  stress  amplitude,  which  exceeds  unity,  is  also  determined  by  the  other  two  stress 
sequences.  Notice  that  as  the  stress  amplitude  in  figure  1 1(a)  grows  without  bound  over  time  for 
a  harmonic  forcing  condition,  the  stress  amplitude  in  figure  11(b)  does  not  exceed  double  the 
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loading  for  a  Heaviside  loading  condition.  Examples  of  such  material  designs  were  given  in  the 
previous  subsection  for  the  four-layer  case. 


(a) 


(b) 


Figure  1 1 .  Stress  time  history  for  a  two-layered  Goupillaud-type  strip  with 

impedance  ratio  ai  =  3,  6\  =  r  =  1,  at  the  middle  of  the  second 
layer  located  at  £  =  3/4  subjected  to  loading  (a)  /(n)  =  sin(nu;)  with 
uj  =  0  =  (b)  f(n )  =  1  for  n  >  0. 


4.3  Natural  Frequencies  of  a  Free-fixed  Non-Goupillaud-type  Layered  Strip  with  Integer 
Layer  Length  Ratios 

The  natural  frequency  results  for  a  free-fixed  Goupillaud-type  layered  strip  with  equal  layer 
lengths  given  in  equations  F-l  through  F-9,  may  be  extended  to  a  free-fixed  non-Goupillaud-type 
layered  strip  with  integer  layer  length  ratios.  We  illustrate  this  with  the  two  examples  included 
below. 

A  three-layered  Goupillaud-type  strip  with  equal  layer  lengths  (a2  =  1  or  ot\  =  1),  represents  a 
two-layered  non-Goupillaud-type  strip  with  1:2  or  2:1  wave  travel  time  ratios  respectively.  The 
case  with  1:2  ratio  requires  that  a2  =  1  and  equation  F-3  becomes 

2  T  2  T  CCi 

cos  —cu  =  —1  or  cos  —uj  =  - - ,  for  a i  >  0.  (65) 

3  3  1 H-  ex  i 

Similar  results  to  equation  65  can  be  obtained  for  the  case  with  2:1  ratio  which  requires  that 

cti  =  1. 

The  frequency  results  for  a  free-fixed  four-layered  Goupillaud-type  strip  with  equal  layer  lengths 
given  in  equation  F-5,  can  be  extended  to  a  free-fixed  three-layered  non-Goupillaud-type  strip 
with  wave  travel  time  ratios  1:2:1  by  choosing  a2  =  1.  With  these  conditions,  equation  F-5 
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becomes 


(1  +  o;i)(l  +  CI3)  cos2  — c 0  —  T3  cos  —  u  +  (r3  —  (1  +  on)(l  +  03)  +  2)  —  0,  (66) 


where  T3  =  (aia3  —  1)  is  chosen  so  that  the  (cosine)  solutions  vary  in  the  interval  [—1, 1].  This 
can  be  achieved  when  T3  =  0,  which  implies  that  (l  +  a1)(l  +  a;3)>2.  Asa  result,  equation  66 
reduces  to: 


0  <  cos2 


(1  +  Qfi)  (1  +  a3)  —  2 
(1  +  ctl)(l  +  0:3) 


(67) 


4.4  Representation  of  the  Stress  Solutions  in  Terms  of  the  Chebyshev  Polynomials 

The  formulas  related  to  the  stress  solutions  in  equations  30,  35,  and  38  consist  of  sums  with  terms 
of  the  form  cos(n0),  sin (n0),  (— l)n,  n  and  (— 1) " (n  +  1).  These  terms  relate  to  the  Chebyshev 
polynomials  of  the  first  and  second  kind  Tn(y )  and  Wn (y) ,  respectively  (see  Schaum  (44)),  as 
shown  below: 


cos(n0)  =  Tn(cos0),  (-ir  =  Tn(-l),  n  =  Wn(l)  -  1,  (-l)n(n  +  1)  =  Wn(-1).  (68) 

In  addition,  from  the  definition  of  the  Chebyshev  polynomials  when  cos  0  7^  0,  one  can  also 
derive  that 


sin(n0)  =  [  IC„(cos0)  —  Tn(cos0)  ]  •  tan0. 


(69) 


As  a  result,  the  stress  formulas  given  by  equations  30,  35,  and  38  can  be  written  in  terms  of 
Chebyshev  polynomials.  We  demonstrate  this  for  the  general  stress  representation  in  equation 
20,  which  can  be  written  in  terms  of  the  Chebyshev  polynomials  as  shown  below: 


Si(n)  = 


for  i  —  1,2,. 


Lt  J 

f(n)  *  [  bi'0Tn(-l)  +  X]  (aiifc  -  6ijfe  tan0fc)  •  Tn(cosdk)  +  bi:kta,n0k  ■  Wn(cosdk)  ]- 

fc= 1 

LfJ 

-/( 0)  •  [  bifiTn(- 1)  +  X]  tan 9k)  ■  Tn(cosdk)  +  bltk  tan0k  ■  Wn( cosdk)  ], 

k= 1 


. .  m  and  n  >  0.  This  representation  is  valid  when  cos  6k  0  0  for  all  1  <  k  < 


(70) 

L?J- 
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5.  Conclusions 


The  resonance  frequency  spectrum  in  equation  15  of  the  discrete  model  was  derived  by 
developing  an  analytical  representation  for  the  stress  solutions  when  a  discrete  sinusoidal  forcing 
function  was  applied  at  one  end  of  the  strip  while  the  other  end  was  held  fixed.  The  approach 
was  based  on  the  ^-transform  methods.  Two  key  elements  of  this  approach  were  that  the 
determinant  of  the  resulting  global  matrix  |  Am|  in  the  /"-space  was  a  palindromic  polynomial  with 
real  coefficients  and  that  its  zeros  were  distinct  and  on  the  unit  circle.  An  important  connection 
was  made  in  equation  15  between  the  angles  (90  =  "  and  {± 9 k}kLi),  representing  the  distinct 
zeros  of  |  Am|  and  the  natural/resonance  frequency  spectrum.  We  showed  that  the  resonance 
frequency  spectrum  in  equation  15  for  the  discrete  model  consists  of  all  the  positive  angles 
coterminal  with  90  =  ^ r,  for  m,  odd  only,  and  {±6^}^ ,  for  any  m.  It  is  important  to  notice  that 
the  positive  angles  coterminal  with  90  =  n  represent  resonance  frequencies  for  any  strip  with  an 
odd  number  of  layers  and  do  not  depend  on  the  material  properties.  According  to  equation  15, 
the  resonance  frequency  spectrum  for  the  discrete  model  depends  only  on  (certain  combinations 
of)  the  impedance  ratios  and  not  on  any  other  parameters  such  as  the  length  of  the  strip,  etc.  The 
design  parameters  representing  such  combinations  for  an  m-layered  strip  were 

m—  1 

Xm- 1  =  n  (1  +  ai)  f°r  2  <  m  <  5  and  Tm_i  for  m  =  4,  5.  This  means  that  as  long  as 

1=1 

m-layered  designs  with  different  impedance  ratios  have  such  design  parameters  in  common,  they 
are  guaranteed  to  have  the  same  resonance  frequency  spectrum  for  the  discrete  model. 

The  natural  frequency  spectrum  of  a  free-fixed  Goupillaud-type  strip  with  up  to  five  layers  was 
first  derived  by  simplifying  and  then  solving  the  frequency  equation.  Independently  of  that,  the 
linear  transformation  in  frequency  in  equation  21  enabled  us  to  heuristically  describe  the  natural 
frequency  spectrum  of  a  free-fixed  m-layered  Goupillaud-type  strip  from  the  discrete  model  in 
equation  15.  Both  approaches  provided  an  easy  way  to  analytically  describe  the  natural 
frequency  spectrum  versus  solving  the  frequency  equation  directly.  The  frequency  spectrum  in 
equations  59  and  appendix  D,  or  appendix  F,  suggests  that  the  natural  frequency  spectrum  of  a 
free-fixed  m- layered  Goupillaud-type  strip  depends  only  on  (combinations  of)  the  impedance 
ratios  and  (it  is  inversely  proportional  to)  the  equal  wave  travel  time  for  each  layer  Designs 
with  these  quantities  in  common  share  the  same  natural  frequency  spectrum  for  the  free-fixed 
boundary  conditions  (see  subsections  3.3.1  and  3.3.2)  and  subsection  4.1  for  illustrations  of  two 
different  ways  of  obtaining  such  designs).  A  possible  design  modification  that  provides  a  desired 
natural  frequency  spectrum  is  also  demonstrated  in  subsection  3.3.4. 
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When  the  resonance  frequency  results  are  combined  with  the  stress  optimization  results  obtained 
in  Velo  and  Gazonas  (40),  we  conclude  that  for  a  given  optimal  design  of  a  Goupillaud-type 
layered  strip  with  one  fixed  end,  depending  on  the  type  of  loading  at  the  other  end,  one  can  get 
either  the  best  or  worst  outcome  when  it  comes  to  the  stress  amplitude.  Although  this  work 
focused  on  a  Goupillaud-type  layered  elastic  strip,  the  frequency  results  were  shown  to  extend  to 
a  non-Goupillaud-type  layered  strip.  The  stress  formulas  involving  trigonometric  sums  were  also 
interpreted  using  the  Chebyshev  polynomials  of  the  first  and  second  kind. 

In  this  work,  we  have  proven  that  the  palindromic  polynomial  with  real  coefficients  |  Am|  has  all 
its  zeros  on  the  unit  circle  for  1  <  m  <  5,  see  appendix  B.  Whether  this  is  true  for  any 
m-layered  Goupillaud-type  strip  or  only  for  certain  classes  of  designs  remains  to  be  investigated. 
Proving  rigorously/theoretically  the  conversion  of  the  resonance  frequency  spectrum  from  the 
discrete  model  in  equation  15  to  the  continuous  model  in  equation  59  also  remains  to  be  shown. 
Finally,  the  experimental  validation  of  the  frequency  results  proposed  here  would  be  a  natural  and 
valuable  extension  of  this  work. 
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Appendix  A.  Stress  Formulas  for  a  Heaviside  Loading 


Some  of  the  definitions  and  formulas  in  appendices  A-E  refer  to  the  results  from  Velo  and 
Gazonas  (40).  These  results  were  obtained  when  studying  stress  propagation  in  a 
Goupillaud-type  layered  strip,  subjected  to  zero  initial  conditions,  a  Heaviside  step  in  stress 
loading  p  =  constant  at  time  t  —  0  at  one  end  and  held  fixed  at  the  other  end.  As  before,  the 

m—  1 

design  parameters  are  given  by  Xm-i  =  II  (1  +  ai)  f°r  2  <  m  <  5,  T3  =  a4a 3  —  1,  and 

i=  1 

r4  =  aia3a4  +  a4a2a4  +  ol\Ol4  +  a2a4  +  a4a3  —  1,  where  Xm-i  >  1  for  m  >  2  and  Tm_!  >  —1 

for  m  =  4,5.  Here  we  recover  from  the  general  formulation  in  equations  5  and  9,  the  stress 
formulas  derived  in  Velo  and  Gazonas  (40)  when  the  strip  is  subjected  to  a  Heaviside  step  in 
stress  loading  p  at  time  t  =  0  on  one  end  and  held  fixed  on  the  other  end.  This  case  is  illustrated 
in  figure  1  for  the  choice  of  f('n)  =  p  for  n  >  0.  Indeed,  in  this  case  we  have 
F(z)  =  Z(f(n ))  =  -p  and  /( 0)  =  p,  which  imply  that  (F(z)  —  /( 0))  =  As  a  result, 

equation  9  becomes 


%m.  — 


|Ai,i| 

|Ai,i| 

Vi z{F(z)  -  /( 0)) 

1  -A-1,2 1 

—  A+2 

7]i  pz 

s 

T— 1 

1 

O 

1 

(~l)1+m\Ahm\ 

(-l)1+m|Alim| 

(A- 1) 


simplifies  to, 


xm(i  = 


PZ  (z—l)\Am\ 


=  pz 


LfJ 


,  Vo_  ,  y-J  ajkz+bjk 

2- 1  ~  2+1  ~  Z-I  22-22COs6»fc  +  l 

k=  1 


(A-2) 


After  applying  the  inverse  ^-transform, 

Si(n)  =  Z~l(xm(i))  =  Z 


L— J 

a-ifiZ  ,  bjfiZ  ai,kz~Jtb’i,kz 

2-1  2+1  2^  22_ 
k= 1 


=  z 


-1 


L— J 

‘i, 0"  _i  bj  Q-  2  \  '  {  Q’i.k'  z(z  COS  OP)  | 

PT  “T"  2+1  l  22-22COS  6k  +  l  “t" 

k=  1  V 


■P 


a 

z— 


2  2  cos  +,  +  I 

bi  /.  ■  2  sin  +. 

22  — 22  COS  0fc  +  l 


•P, 


(A- 3) 
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the  stress  formulas  follow, 


Si(n) 


LfJ 

ciifl  +  bifi(-l)n  +  ai,k  cos  (nOk)  +  bi>k  sin  (nOk) 
k= 1 


■P, 


for  the  choices  of  the  coefficients  given  below 


a*  k  cos  6k  +  b*  k  rn 

ai,k  =  aitk,  bt.k  =  — i for  k  =  1, . . . ,  |  — J . 

sm  u  i ^ 


(A-4) 


(A- 5) 


As  before,  sin  6k  ^  0  for  all  k  —  1, . . . ,  |_y  J  and  0  =  0  for  m  even.  In  addition,  the 
trigonometric  sums  in  the  stress  equations  A-4,  allow  us  to  express  them  in  terms  of  the 
Chebyshev  polynomials  of  the  first  and  second  kind  Tn(y )  and  Wn(y),  by  applying  equation  68 
and  equation  69  as  follows 


Si(n) 


LfJ 

lift  +  bifiTn(-l)  +  Y,  ~  bi,k  tan  8k)  •  T„(cos  8k)  +  bl>k  tan  0k  ■  Wn(cos6k) 
k— l 


■P, 


(A-6) 


for  i  —  1,  2, . . .  m  and  n  >  0.  This  representation  is  valid  when  cos  9k  ^  0  for  all  1  <  k  <  [yj . 
The  coefficients  ahj,  bij,  as  well  as  the  angles  6k  depend  only  on  the  impedance  ratios  a,, 
i  =  1,  2, . . . ,  m,  j  =  0, 1, . . . ,  [f  J  and  k  =  1, . . . ,  [f  J . 
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Appendix  B.  Zeros  of  the  Determinant  |  Am 


The  expressions  for  the  determinant  |  Am|  using  the  recursive  equation  10  for  1  <  m  <  5  are: 

| Ai|  —  z  +  1,  |A2|  =  (z  +  l)2  -  r}2rftaiz  =  z2  -  2(x^~2)^  +  1, 

\As\  =  (z  +  l)-[z2-%^z  +  l],  \A4\=z4-^z3  +  ^f±^z2-^z  +  l,  (B-l) 

|  A5  |  =  (Z  +  1)  ■  [  ^  _  4T4.3  +  2(4T4-X4+8)  2  „  _|_  \  ]. 

1  1  V  ’  L  X4  X4  X4  1 

The  determinant  |  Am|  is  a  palindromic  polynomial  with  real  coefficients,  as  expected.  All  the 
zeros  of  |  Am|  for  1  <  m  <  5  are  on  the  unit  circle  and  distinct  based  on  the  following  criteria. 


1.  The  zeros  of  a  quadratic  factor  with  real  coefficients  of  the  form  z2  +  fiz  +  1  are  distinct 
(complex  conjugate)  on  the  unit  circle  iff  yu2  —  4  <  0  or  equivalently  iff  |/i  <  2.  When  m  =  2,3 

<  1,  which  is  true  because  Xm-i  >  1. 


we  have  that 


<  2 


Xm  —  1  ~2 

Xm  —  1 


2.  The  zeros  of  a  quartic  factor  with  real  coefficients  of  the  form  z4  +  /j,z3  +  uz2  +  /iz  +  1  are 


distinct  (complex  conjugate)  on  the  unit  circle  2  iff  [  fi2  —  Av  +  8  >  0  and 


±\/ /i2— 4^+8  | 


<  2 


When  m  —  4, 5  we  have  the  condition: 


f!2  -  Av  +  8  >  0  <^>  [  (Xm- 1  -  rm_i)2  -  4xm-i  ]  =  [  (Tm-1  -  2)2  -  4(1  +  Tm_i)  ]  >  0,  where 
"f'm—  1  Xm—1  fm-1’ 


The  inequality  follows  from  the  following  relations: 


(T 3  —  2)2  —  4(1  +  T3)  >  («i  +  a3)2  —  4aia3  =  (o>i  —  a3)2  >  0  for  m  =  4,  and 
(T4  —  2)2  —  4(1  +  T 4)  >  {oil  "T  ot2  +  0:40:2  +  o:3  +  0:4  +  oi3oi 4)“  —  4(1  +  ^)  > 
(—0:4  —  0:2  —  0:40:2  +  a3  +  0:4  +  o:3o:4)2  >  0  for  m  =  5. 


The  condition 


~n±  4t2— 4^+8  2  r">-l±V  (Xm-l-rm_i)2-4xm-l 

2  Xm  —  1 

Xm- 1  >  1  >  0  and  rm_!  +  1  >  0  for  m  =  4,  5. 


<  1,  which  is  true  because 


^This  criterion  is  derived  after  dividing  the  quartic  polynomial  by  z 2 

l 


and  then  applying  the  transformation  y  = 
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Appendix  C.  Definition  of  the  Matrix  Dm 


The  following  equation  is  the  definition  of  matrix  Dm: 


Dr 


F 

C2 

0 


0 

0 


Bi  0  ... 

F  B2  0 
C3  F  B3 


0  Ci  F  Bi  0 


0 

0 

0 

0 

0 


0  Cm_i  F  Bm_i 
0  Cm  F 


(C-l) 


-  2mx2m 


where  F  =  I  +  RT,  B;  =  ^^1,  Q  =  Cm  =  4RT,  1  <  i  <  m  -  1. 

Here  8  represents  any  of  the  angles  9j,  j  —  0, . . . ,  |_yj ,  as  applicable;  I  represents  the  2  x  2 


identity  matrix;  and  RT  represents  the  transpose  of  the  rotation  matrix  R  = 


cos  9  —  sin  9 
sin  9  cos  9 
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Appendix  D.  Necessary  Relations  Among  the  Base  Angles  and  Design 
Parameters 


The  angles  identified  here  for  a  strip  with  up  to  five  layers  relate  to  the  zeros  of  |  Am|,  already 
shown  to  be  on  the  unit  circle  in  appendix  B.  The  2  \Jf\  order  equation  for  finding  the  non-trivial 
zeros  of  |  Am|  in  equation  13  is  reduced  to  an  [f\  order  equation  for  the  cosine  of  the  base 
angles.  The  inequalities  satisfied  by  the  design  parameters  Xm-i  for  m  =  2,3  and  Xm-i,  rm_i 
for  m  =  4,  5  in  appendix  B,  validate  the  existence  of  these  angles. 


The  two-layer  case.  The  condition  |D2 


0  equivalent  to 


cos  9  — 


XI -2 

Xl 


2 


0,  implies  that 


T  <  cos  9  =  — — -  <  1 
Xi 


0  <  9i  =  cos 


-l  ( Xi  ~  2 
Xi 


<  7T, 


(D-l) 


The  three-layer  case.  The  condition  |D3|  =  0  is  equivalent  to 


cos  9  =  — 1  or  (  cos  9  —  )  =  0 


,  which  based  on  the  fact  that 
X2  =  (1  +  «i)(l  +  ol 2)  >  1,  leads  to  the  relations 


Vo  —  2 

cos#  =  — 1  or  —  1  <  cos 9  =  — -  <  1, 


and  the  corresponding  angle  solutions 

90  =  cos_1(— 1)  —  it,  0  <  9i  —  cos 


X2 


_i  1  X  2  —  2 


(P-2) 


X2 


<  7T. 


(D-3) 


The  four-layer  case.  The  condition  |D4|  =  0  «=>  [x3  cos2  9  -  2T3  cos  9  +  (2T3  -  *3  +  4 )]2  =  0, 

leads  to  the  quadratic  equation: 


2  n  2T3 

cos  9 - cos  9  + 

X3 


2r3  -  xz  +  4 

Xs 


=  0, 


(D-4) 


and  implies  the  following  base  angle  solutions  for  any  given  set  of  the  design  parameters  %3  and 

r3, 


0  <  9\  2  =  cos 


-1 


r3  ±  y/(x3  -  rs) 2  -  4x3 


<  7T. 


X3 


(P-5) 
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The  five-layer  case. 

|D5|  =  0  <^>  [  cos  6  =  — 1  or  [  \4  cos2  6  —  2r4  cos  6*  +  (2r4  —  \4  +  4)  ]2  =  0  ],  leads  to  the 
equations: 

2r4  2r4  —  y4  -t-  4 

cos  9  =  —l  or  cos2  6 - cos 9  H - - - =  0, 

X4  X4 


and  the  corresponding  angle  solutions: 


r4)2  -  4\4  j 


for  any  given  set  of  the  design  parameters  X4  and  h4. 


9 0  =  COS  i(  — 1)  =  7T,  0  <  9 1)2  = 


=  cos 


-1 


r4  ±  y/(x 4  - 


X4 


(D-6) 


(D-7) 
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Appendix  E.  Optimal  Base  Angles 


For  an  m- layered  Goupillaud-type  elastic  strip  subjected  to  a  Heaviside  loading  at  one  end  and 
held  fixed  at  the  other  end,  the  optimal  designs  have  the  smallest  stress  amplitude  of  double  the 
loading.  These  optimal  designs  are  characterized  by  the  following  optimal  values  of  their  base 
angles  for  j  =  2,3,4,...: 


Table  E- 1 .  Optimal  base  angle  values  for  the  m-layer  case. 


m 

2 

3 

4 

5 

Optimal  Base  Angles 

II 

a  7 r 

171, opt  —  2j—l 

0  77 
°l,opt  —  2j  ’ 

@2, opt  —  'K  2j 

a  _  7r 

(71, opt  —  2j+l  ! 
a  _  (2j  — 1)tt 

(72, opt  —  2.7+1 
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Appendix  F.  Solutions  of  the  Frequency  Equations  for  a  Free-fixed 
m-layered  Goupillaud-type  Elastic  Strip  (2  <  m  <  5) 


F-l  Natural  Frequency  Spectrum  for  the  Free-fixed  Two-layer  Case 


After  a  few  mathematical  manipulations  involving  trigonometric  identities,  the  frequency 
equation  57  for  m  =  2  becomes 

Ot\  —  1  Vy  —  2 

costcu  =  - - =  - - ,  where  Xi  —  (1  +  «i)  >  1  or  ay  >  0.  (F-l) 

ai  +  1  xi 


This  suggests  the  following  natural  frequency  spectrum 


Vh 


1 

r 


+  2Zi7t 


1 

r 


+  2(Z2  +  l)71" 


(F-2) 


for  Zi ,  Z2  =  0,1,2,....  Here  a  i  represents  the  impedance  ratio  between  the  first  and  the  second 
layer,  while  Xi  —  1  +  «i-  When  both  layers  are  occupied  by  the  same  material  (o:i  =  1), 
equation  F-l  recovers  a  known  literature  result  for  the  natural  frequency  spectrum  of  a 
homogeneous  strip  of  length  r  given  in  Graff  (33).  In  this  case  the  frequency  equation  F-l 
becomes  cos  tuj  =  0,  which  implies  the  natural  frequency  spectrum  u>i  =  J,_  lj  -|,/  =  0,1,2.... 
Notice  that  equation  F-l  implies  that  cos  tuj  ^  ±1. 


F-2  Natural  Frequency  Spectrum  for  the  Free-fixed  Three-layer  Case 

Using  the  fact  that  |  Am|  =  |  Am|,  the  frequency  equation  57  for  m  =  3  becomes  equivalent  to 
solving  the  following  equations: 

2  T  2  T  X2  —  2 

cos  — lo  =  —1  or  cos  — u  =  - - ,  where  %2  =  (1  +  «i)(l  +  a2)  >  1.  (F-3) 

3  3  X2 

From  equation  F-3,  the  frequency  spectrum  for  the  three-layer  is 


~  (2^o  +  1)^, 


(F-4) 


= 


2  T 


COS 


-1  1  X2-2 
X2 


+  2li1T 


,  = 


2  T 


COS 


-1  l  X2-2 
X2 


+  2  (/2  +  1)7T 
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where  lQ,  h,  Z2  =  0,1,2,.... 


F-3  Natural  Frequency  Spectrum  for  the  Free-fixed  Four-layer  Case 


Use  of  the  double  angle  formula  for  the  cosine  for  the  four-layer  case  (m  —  4),  results  in 


|A4|  =  \  •  [x3  COS2  \00 


2r3  cos  \oo  +  (2r3  —  X3  +  4)]  =  0,  or  equivalently 


r  2r, 


r  2r3-Xs  +  4 


cos  -o o - COS  —00  + 

2  X3  2  X3 


=  0. 


(F-5) 


For  any  given  values  of  the  design  parameters  X3  and  T3,  the  frequency  spectrum  is: 

1  ( r3±V(X3-r3)2-4x3\  +  2lJ  , 


Uh  = 


Ul2  = 


cos 


X3 


h,  h  —  0, 1, 2  ... , 


(F-6) 


cos” !  ^r3±V(X3-r3)2-4x3^  +  2(Za  +  1)7r 


When  the  first  two  layers  are  occupied  by  one  material  (an  =  1)  and  the  other  two  are  occupied 
by  another  material  {a.2  1 ,  a3  =  1),  the  four-layer  case  reduces  to  the  two-layer  case  and 

equation  F-5  reduces  to  equation  F-l.  From  equation  F-5  if  follows  that 


9  T  Ol2 

COS"  —00  =  - . 

2  a2  +  1 


(F-7) 


which  after  applying  the  double  angle  formula  for  cosine  becomes  consistent  with  equation  F-l, 

a2  -  1 


COS  TOO  = 


OL2  +  1' 


(F-8) 


F-4  Natural  Frequency  Spectrum  for  the  Free-fixed  Five-layer  Case 

Similar  to  the  previous  cases,  using  symbolic  algebra  software  and  trigonometric  identities  we 
obtain  for  rn  —  5  that  |  A5|  =  cos  ^  •  [y4  cos4  \oo  —  (x4  +  T4)  cos2  +  (r4  +  1)]  =  0.  After 
using  the  double  angle  trigonometric  formulas  we  have  that 


2  r 


21- 


cos  —oo  =  —  1  or  cos  —oo  — 


2r4 


2  r 


_  -  - _  -  cos  —oo  + 

5  5  xr  5  X4 

The  frequency  spectrum  can  then  be  obtained  for  any  given  values  of  x4  and  F4. 


2T4  —  X4  +  4 


=  0. 


(F-9) 
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Appendix  G.  Material  Properties 


The  elastic  modulus  E  and  density  p  of  the  materials  used  in  subsection  4.1  (see  Ashby  (45))  are 
as  follows: 

•  Aluminum  alloy:  E  —  70  GPa  and  p  =  2,  500  kg/m3, 

•  Brass:  E  —  20  GPa  and  p  =  984.19  kg/m3, 

•  Bronze:  E  —  50  GPa  and  p  =  7, 189.35  kg/m3, 

•  Copper  alloy:  E  —  10  GPa  and  p  =  515.19  kg/m3, 

•  Glass  fiber  reinforced  material:  E  =  30  GPa  and  p  =  1, 130  kg/m3, 

•  Lead:  E  =  14  GPa  and  p  =  11,  340  kg/m3, 

•  Nickel  alloy:  E  =  200  GPa  and  p  =  348.18  kg/m3, 

•  Tungsten  alloy:  E  =  275  GPa  and  p  =  19,  610  kg/m3, 

•  Steel:  E  =  200  GPa  and  p  =  8,  000  kg/m3. 
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